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ABSTRACT 

This  report  describes  the  structure  and  use  of  Monte  Carlo  programs  that  calculate 
the  transport  of  proton  beams  through  extended  media.  Although  more  generally 
explicable,  the  programs  have  been  designed  to  deal  with  the  penetration  of  50-  to 
250-MeV  beams  through  water  phantoms.  The  Monte  Carlo  model  takes  into 
account  multiple-scattering  deflections  and  energy-loss  straggling  due  to  Coulomb 
interactions  of  protons  with  atoms  and  orbital  electrons.  Nonelastic  nuclear 
interactions  are  treated  as  an  absorptive  effect.  The  PTRAN  system  at  present 
consists  of  several  cross-section  preparation  programs  and  two  main  codes, 
PTRAN3D  and  PTRANID.  PTRAN3D  ^plies  to  an  incident  narrow  pencil 
beam,  and  calculates  (a)  the  deposition  of  energy  as  function  of  depth  and  radial 
distance  from  the  beam  axis,  and  (b)  and  the  energy  spectra  of  the  primary 
protons  as  function  of  depth.  Program  PTRANID  is  a simplified  version  which 
runs  faster  and  omits  the  calculation  of  the  radial  distribution  of  energy  deposition. 


*Contractor,  work  done  under  NIST  contract  50SBN2C7042. 


1.  Introduction 


This  report  describes  a system  of  Monte  Carlo  programs  for  calculating  the  penetration, 
diffusion  and  slowing  down  of  proton  beams  in  an  extended,  homogeneous  medium.  The 
objeaive  is  to  provide  information  useful  for  treatment  planning  and  dosimetry  in  proton  therapy. 
Attention  is  therefore  focussed  on  the  penetration  through  water  of  proton  beams  with  initial 
energies  from  50  MeV  to  250  MeV. 

The  Monte  Carlo  model  is  based  on  the  condensed-random-walk  method  (Berger,  1963), 
and  takes  into  account  the  following  types  of  events  occurring  in  successive  short  track  segments: 
(a)  energy-loss  straggling  in  Coulomb  collisions  with  atomic  electrons,  (b)  multiple-scattering 
deflections  due  to  elastic  scattering  by  atoms,  and  (c)  energy  losses  in  nonelastic  nuclear 
reactions.  The  Monte  Carlo  model  can  be  applied  to  a variety  of  transport  problems.  In  this 
report,  programs  are  described  which  calculate  the  spatial  distribution  of  energy  deposition  (in 
one  or  three  dimensions),  as  well  as  energy  spectra  of  protons  at  various  depths  in  a phantom. 

In  Section  2,  procedures  for  the  simulation  of  proton  tracks  are  described.  In  Section  3, 
data  preparation  programs  are  described  which  facilitate  the  sampling  of  energy  losses  from  the 
distribution  of  Vavilov  (1957)  and  the  sampling  of  angular  deflections  from  the  distribution  of 
Moli^re  (1948).  With  these  procedures  one  can  construct  Monte  Carlo  programs  for  treating  a 
variety  of  proton  transport  problems.  Two  such  programs  are  described  in  this  report.  Section  4 
deals  with  program  PTRAN3D,  which  calculates  the  deposition  of  energy  as  function  of  depth 
and  of  the  radial  distance  from  a narrow  pencil  beam,  and  also  proton  spectra  as  functions  of 
depth.  Section  5 describes  a simpler  one-dimensional  program  PTRANID,  which  omits  the 
calculation  of  the  radial  energy-deposition  distribution.  Section  6 discusses  the  random-number 
generators  used  in  the  PTRAN  programs.  Section  7 provides  information  about  the  dependence 
of  the  results  on  the  step-size  of  the  condensed-random  walk  model.  In  Section  8,  the  program 
and  data  files  are  listed  which  are  included  in  the  PTRAN  system  at  the  present  time. 

For  the  purpose  of  calculating  the  spatial  distribution  of  energy  deposition,  the  transfer  of 
energy  from  the  proton  beam  to  the  medium  can  be  treated  as  a two-stage  process.  In  the 
important  first  stage  one  considers  the  spatial  distribution  of  (1)  proton  energy  losses  that  occur 
as  the  result  of  Coulomb  interactions  with  atoms  or  molecules,  and  (2)  of  the  removal  of  protons 
from  the  beam  due  to  nonelastic  nuclear  interactions.  In  the  second  stage  one  takes  into  account 
the  further  spatial  transport  of  energy  by  secondary  electrons  from  ionization  events,  and  by 
charged  particles,  neutrons,  or  gamma  rays  from  nuclear  reactions.  The  Monte  Carlo  programs 
described  in  this  report  deal  only  with  the  first  stage. 

The  ranges  of  the  secondary  electrons  are  exceedingly  short  compared  to  the  ranges  of  the 
primary  protons.  Therefore  energy  transport  by  secondary  electrons  has  an  effect  on  depth  dose 
curves  only  at  very  shallow  depths,  where  it  gives  rise  to  a rapid  dose  buildup.  In  regard  to 
nonelastic  nuclear  reactions,  one  must  estimate  the  fraction  of  the  energy  lost  by  protons  that 
esc^e  in  the  form  of  neutrons,  x-rays,  or  fast  protons,  and  the  fraction  that  can  be  considered  to 
be  absorbed  locally.  It  is  also  of  interest  to  calculate  the  energy  degradation  spectra  of  the 
secondary  heavy  charged  particles.  Work  on  these  topics  is  in  progress  at  NIST  by 
S.  M,  Seltzer,  and  the  results  will  be  integrated  into  the  analysis  of  output  from  the  PTRAN 
codes. 
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2.  Monte  Carlo  Model 


2.1,  Schematization. 

For  the  purpose  of  simulation,  proton  tracks  are  divided  into  many  short  segments. 

These  segments  are  also  called  the  "steps"  of  a condensed  random  walk.  Energy-losses  in 
successive  steps  are  sampled  from  the  Vavilov  distribution  (Vavilov,  1957).  The 
multiple-scattering  angular  deflection  in  successive  steps  are  sampled  from  the  Molifere 
distribution  (Moli6re,  1948).  Nonelastic  nuclear  reactions  are  treated  an  absorptive  effect. 
Appendices  1 and  2 list  the  formulas  used  for  calculating  the  Molibre  and  Vavilov  distributions. 

In  condensed-random-walk  model  used  in  the  present  work,  a proton  track  (Monte  Carlo 
history)  is  described  by  the  following  array: 


So» 

Si. 

S2,  — , 

Sn» 

Eo» 

El, 

E2, 

En. 

Xl. 

X2, 

Xn. 

yo» 

yi. 

y2.  ••••. 

yn» 

Zi, 

Z2,  

Zn. 

^0’ 

^1. 

$2,  — , 

^1, 

<P2.  •••'. 

^n» 

Wo, 

Wi, 

W2,  ...., 

w 

(2.1) 


The  index  n pertains  to  the  track  characteristics  at  the  end  of  the  n^  step,  s^  is  the  path 
length  traveled  by  the  proton  since  its  entry  into  the  medium,  and  As^  = s^^.^  — s^  is  the  size  of 
the  n^  step.  E„  is  the  kinetic  energy;  x„,  y^  and  z^  are  spatial  coordinates;  6^  and  v?n  specify  the 
direction  of  motion  in  a spherical  coordinate  system  whose  polar  axis  is  the  z-axis;  is  a 
survival  weight  factor  which  represents  the  probability  that  the  proton  has  not  been  absorbed,  in 
path  length  Sn,  by  a nonelastic  nuclear  reaction.  The  index  o pertains  to  the  initial  conditions,  at 
the  p)oint  of  entry  of  the  proton  into  the  medium,  so  that  Sq  = 0 and  Wq  = 1. 

2.2.  Choice  of  Step  Sizes. 

Several  considerations  enter  into  the  choice  of  the  step  sizes  of  the  condensed  random 
walk.  In  order  to  reduce  the  computing  time,  the  steps  should  be  as  large  as  possible.  In  order 
to  reduce  the  error  resulting  from  the  neglect  of  individual  collisions,  the  steps  should  be  as  short 
as  possible.  The  steps  must  be  sufficiently  long  so  that  the  Vavilov  energy-loss  distribution  and 
Moli^re  multiple-scattering  distribution  are  both  applicable.  On  the  other  hand,  the  steps  should 
be  short  enough  so  that  the  energy  loss  per  step  is  a small  fraaion  of  the  proton  energy  at  the 
beginning  of  the  step. 

In  the  implementation  of  the  condensed  random  walk  model  in  the  PTRAN  system,  the 
choice  of  step  sizes  is  made  via  the  adoption  of  an  energy  grid  at  which  the  various  cross  sections 
and  probability  distributions  are  evaluated.  This  arrangement  will  be  discussed  in  Section  3.1, 
Some  numerical  experiments  were  carried  out  to  explore  the  effects  of  changing  the  step  sizes, 
and  will  be  discussed  in  Section  7. 
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2.3.  Energy  Loss  and  Change  of  Direction. 

The  energy  loss  in  the  n***  step,  — E„+j,  is  sampled  from  the  Vavilov  distribution. 
The  multiple-scattering  angular  deflection  in  the  n***  step  is  specified  in  terms  of  two  angles, 
d'  and  <p\  in  a coordinate  system  whose  polar  axis  coincides  with  the  direction  of  motion  of  the 
proton  at  the  beginning  of  the  segment.  The  angle  d'  is  sampled  from  the  Moli^re  distribution, 
and  the  azimuthal  angle  from  a uniform  distribution  between  —180  and  180  degrees.  The 
direction  cosines  of  the  proton  track  at  the  end  of  the  n***  step  are  given  by 


sm9„.,  cos(0„.i 

sind'  cosv?' 

= R • 

sind'  sin^' 

cos9„.i 

cos^' 

» t 

where  R is  the  rotation  matrix 

COS0„  COS<PJ^ 

Sin^n  COSV?n 

R 

= 

cos^n  sin^n 

COS^n 

sin^n  sin<«>n 

-sin9„ 

0 

COS^n 

In  the  PTRAN  codes,  Eqs  (2.2)  and  (2.3)  are  replaced  by  the  equivalent  equations 


cos0„^i 
sin((p„^i  - (Pn) 

- ^n) 


= cos^pcos^'  - sin^n  sin0' cos^' 
sin0'sin^' 

cosd'  - COS0„  COS^n+i 
sin0„sin0„,i 
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(2.2) 


(2.3) 


(2.4) 

(2.5) 

(2.6) 


When  sin^n  sin^^+i  = 0,  then  ^^+1  - <Pn  = 

2.4.  Change  of  Position. 

Let  Ax',  Ay',  and  Az'  denote  the  spatial  displacements  or  the  proton  in  the  n^  step,  in  a 
Cartesian  coordinate  system  whose  axis  coincides  with  the  direction  of  motion  at  the  beginning  of 
the  step.  Taking  into  account  multiple  scattering,  these  displacements  are  calculated  from  the 
following  ^proximate  expressions: 


3 


Ax' 

A>'' 

Az' 


1^" 


sin0'  cos^'  + kj 
sin^'  sin^p'  + ky  (?/6)^^ 


1 ASn  (1  + COS0')  . 


(2.7) 

(2.8) 

(2.9) 


Here  6'  and  <p'  are  the  multiple-scattering  angular  deflections,  and  kj^  and  ky  are  random  variables 
distributed  according  to  a Gaussian  distribution  with  zero  mean  and  unit  standard  deviation.  The 

quantity  ^ is  the  mean-square  multiple-scattering  deflection  which  is  approximated  by 
2 2 

Xc(B-1.2),  where  Xc  B are  variables  in  Moli6re’s  multiple-scattering  theory  (see 
Appendix  1). 

The  corresponding  spatial  displacements  in  a fixed  coordinate  system  are  related  to  Ax', 
Ay',  and  Az'  by  the  rotation  matrix  R from  Eq  (2.3): 


^n+l 

Ax' 

yn-.!  - yn 

= R • 

Ay' 

Zn.l  - Zn 

Az' 

(2.10) 


2.5.  Nuclear  Survival  Weights. 

The  weight  factor  is  treated  in  PTRAN  on  the  basis  of  the  continuous-slowing-down 
^proximation,  assuming  that 


W„(T)  = exp 


dT' 

S(T') 


(2.11) 


where  is  a nuclear  absorption  coefficient  that  represents  the  probability,  per  unit  path  length, 
of  a nonelastic  nuclear  reaction.  This  coefficient  is  given  by 


N 


M„uc(T)  = a(T)  , 


nuc 


(2.12) 


where  is  Avogadro’s  number,  A the  molecular  weight,  and  a the  total  non-elastic  nuclear 
reaction  cross  section. 
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3.  Data  Preparation  Programs 


Proton  stopping  powers  and  ranges  were  calculated  with  the  program  PSTAR  (Berger, 
1992)  which  generates  results  consistent  with  those  tabulated  in  Report  49  of  the  International 
Commission  on  Radiation  Units  and  Measurements  (ICRU,  1992).  Stopping  powers  and  ranges 
in  water  are  given  in  table  1 . 

The  nonelastic  nuclear  cross  section  for  hydrogen  is  negligible.  The  adopted  cross 
section  for  oxygen,  shown  in  figure  1,  is  based  on  a fit  to  experimental  results  of  Renberg  et  al. 
(1972)  and  Carlson  et  al.  (1975).  In  the  energy  region  from  threshold  up  to  10  MeV  the 
energy-dependence  of  the  adopted  cross  section  is  modeled  after  theoretical  results  (S.  M.  Seltzer, 
private  communication)  obtained  with  the  GNASH  code  of  Young  et  al.  (1990). 

3.1.  Program  PARAM4. 

A sequence  of  decreasing  energies  Tq,  Tj,  T2,  •••,  T^,, ...,  Tj^  is  chosen  such  that  the 
difference  T^,  — has  either  a constant  value  AT,  or  the  value  kT^,  whichever  is  smaller. 
In  the  exploratory  calculations  described  in  this  report,  five  such  energy  grids  were  used,  with 
AT  = 4,  2,  1,  0.5  and  0.25  MeV,  and  with  k = 0.05  in  all  cases. 

For  each  energy  interval  (T^,  T^+i),  a path  length  6s^,  is  calculated  in  the  continuous- 
slowing-down  ^proximation: 


where  S(T)  is  the  stopping  power.  The  set  of  path  lengths  Ss^,,  m = 1, M,  will  be  referred  to 
as  a step-size  grid.  Actual  step-sizes  are  chosen  in  the  PTRAN  programs  by  interpolation  in  a 
table  of  6s„  vs.  T^. 

Various  parameters  of  the  Molibre  and  Vavilov  distribution  are  calculated  in  PARAM4 
for  path  lengths  6s^,  using  the  formulas  in  Appendices  1 and  2.  The  Vavilov  parameters  are 
evaluated  for  an  initial  energy  T^,.  The  Molifere  parameters  are  evaluated  at  an  intermediate 
energy  (Tn,  + T„+i)/2,  in  order  to  account  for  the  energy  loss  along  the  track.  PARAM4  also 
tabulates  the  values  of  the  proton  stopping  power  and  range,  the  nuclear  attenuation  coefficient, 
and  the  survival  weight  factor  W,  at  energy  T^^. 

Input  Parameters.  The  user  is  prompted  to  provide,  from  the  keyboard,  the  grid 
parameters  AT  and  k,  and  the  starting  energy  Tq  and  length  M of  the  grid.  On  the  monitor 
screen  it  is  indicated  how  these  quantities  can  be  chosen  to  obtain  the  five  grids  mentioned  above. 
The  user  is  free,  however,  to  m^e  different  choices.  The  user  must  also  specify  the  names  of 
the  ouftiut  files,  and  suggested  names  are  supplied. 
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Input  Data  Files.  The  following  files  are  required: 

COMPOS. WAT  Composition  data  for  water 

STOPRANG.WAT  Stopping  power  and  range  table  for  water 

OXFIT.DAT  Cross  section  for  non-elastic  reactions  of  protons  with 

oxygen  nuclei 

PARAM4  could  be  adapted  to  materials  other  than  water.  Small  modifications  of  the  source  code 
would  be  necessary,  which  are  indicated  in  the  program  listing. 

Output  from  PARAM4:  Three  ouq)ut  files  are  generated: 

1)  Ouqiut  table  for  inspection  by  the  user,  with  a suggested  name 
PARAM4.TBj,  where  j identifies  the  step-size  grid. 

2)  A file  of  data  to  be  used  as  input  for  programs  VPREP4  and 
MPREP4,  with  the  suggested  name  PARAM4.ARj. 

3)  A file  of  data  to  be  used  as  input  to  the  main  Monte  Carlo  program 
(such  as  PTRANID  or  PTRAN3D),  with  the  suggested  name 
PARAM4.PTj. 

Output  table  PARAM4.TBj  lists  v^ious  quantities  pertaining  to  the  Molifere  and  Vavilov 
distributions,  as  well  as  nuclear  attenuation  coefficients  and  survival  factors.  By  inspecting  this 
table,  the  user  can  determine,  for  example,  whether  the  values  of  the  Moli6re  parameter  B,  as 
required,  are  greater  4.5;  or  whether  the  values  of  the  Vavilov  parameter  are,  as  required,  much 
larger  than  the  mean  excitation  energy  of  the  medium.  An  excerpt  from  output  file 
PARAM4.TB4,  for  Grid  4,  is  shown  in  table  2. 

3.2.  Program  VPREP4  and  Related  Programs. 

Arrays  for  Alias  Sampling  Method.  In  PTRAN,  energy  losses  are  sampled  from  the 
Vavilov  distribution  by  a two-step  procedure.  First,  a random  selection  is  made  of  a bin  in  a 
histogram  of  the  Vavilov  distribution,  and  then  a value  of  the  scaled  energy  loss  is  chosen  at 
random  within  this  bin. 

VPREP4  calculates  the  Vavilov  distribution  according  to  Eq  (A2.6)  as  a 

function  of  the  scaled  energy-loss  variable  X,  the  proton  velocity  and  the  skewness  parameter  k 
(see  Appendix  2).  The  integral  in  Eq  (A2.6)  is  evaluated  by  an  adaptive  numerical  quadrature 
code  (Kahaner  et  al.,  1989).  The  distribution  is  then  converted  into  a histogram  with  a constant 
bin  width  AX  = 0.1. 

The  random  selection  of  a histogram  bin  is  done  by  the  alias  sampling  method  (Walker, 
1974;  Kronmal  and  Peterson,  1979),  which  has  the  advantage  that  a single  random  number  and  a 
single  comparison  are  sufficient.  VPREP4  calculates,  from  the  Vavilov  histogram,  two  auxiliary 
arrays  which  are  required  by  the  alias  sampling  method.  This  is  done  with  subroutine  ALIAS 
which  is  an  adaptation  of  a program  developed  by  Kronmal  and  Peterson  (private 
communication). 
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Output  Options.  There  are  two  options,  one  for  production  and  the  other  one  for  testing: 

1)  Output  for  use  in  PTRAN  programs.  The  arrays  needed  for  the  alias  sampling  method 
are  calculated  for  many  grid  intervals  between  specified  lower  and  upper  limits.  The  ou^ut  is 
stored  in  a large  formatted  file  with  the  suggested  name  VPREP4.ARj,  where  j indicates  the  grid. 
In  order  to  reduce  the  time  needed  to  enter  these  data  into  the  PTRAN  programs,  a conversion 
program  VCON4  is  used  to  generate  from  VPREP4.ARj  a corresponding  binary  file  with  the 
suggested  name  VPREP4.BRj.  Such  a file  can  provide  input  for  Monte  Carlo  calculations 
involving  proton  beams  with  many  different  energies. 

2)  Output  for  Testing.  The  output  file  contains  the  two  alias-method  arrays,  and  the 
underlying  probability  histogram,  for  a single  grid  interval.  These  data  can  be  used  as  input  for 
an  auxiliary  program  called  VSAMP4,  which  samples  energy  losses  by  the  alias  method.  The 
sampled  histogram  can  be  compared  with  probability  histogram  to  verify  that  the  sampling 
procedure  works  correctly. 

With  this  option  it  is  also  possible  to  generate  an  output  file  which  contains  the  Vavilov 
distribution  (rather  than  a histogram)  as  a function  of  the  scaled  energy-loss  parameter  X. 

Vavilov  distributions  for  grid  intervals  with  starting  energies  T^  = 25,  50,  100  and  200  MeV 
(for  Grid  4)  are  illustrated  in  figure  2.  Program  VSUM4  uses  numerical  quadrature  to  obtain  the 
mean  value  and  variance  of  the  Vavilov  distribution.  This  allows  the  user  to  check  that  the  mean 
value  equals  T^^  - and  that  the  variance  equals  the  theoretical  variance  from  Eq  (A2.13). 

Finally,  with  option  2 it  is  also  possible  to  omit  the  Blunck-Leisegang  correction,  setting  e = 0 
in  Eq  (A2.8)  of  Appendix  2.  For  the  grids  used  in  this  paper,  this  omission  would  have  very 
little  effect  on  the  calculated  Vavilov  distributions. 

3.3.  Program  MPREP4  and  Related  Programs. 

MPREP4  calculates  the  Moli&re  distribution  for  each  path  length  6Sn,,  using  the  equations 
in  Appendix  1.  The  distribution  is  evaluated  at  156  values  (between  0 and  40)  of  the  reduced 
Molitre  angle  t?.  The  spacing  of  the  angles  is  non-uniform,  and  increases  with  increasing  t?. 

The  set  of  t?- values  is  contained  in  file  THSET2.  A probability  histogram  with  155  non-uniform 
bins  is  then  calculated,  and  is  used  to  obtain  the  two  auxiliary  arrays  for  die  alias  sampling 
method.  There  are  two  output  options: 

1)  Output  for  use  in  PTRAN  programs.  The  output  file  contains  the  two  arrays  for  the 
alias  sampling  method,  for  the  entire  grid.  The  output  file  is  a large  formatted  file  with  the 
suggested  name  MPREP4.ARj.  Program  MCON4  is  used  to  generate  a corresponding 
unformatted  (binary)  file  MPREP4.BRj. 

2)  Output  for  Testing.  The  output  includes  the  two  alias  method  arrays,  and  the 
underlying  probability  histogram,  for  a single  grid  interval.  These  data  can  be  used  as  input  for 
an  auxiliary  program  MSAMP4,  which  samples  angular  deflections  by  the  alias  method.  The 
sampled  histogram  can  be  compared  with  input  probability  histogram,  to  verify  that  the  sampling 
procedure  works  correctly. 
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4.  Monte  Carlo  Program  PTRAN3D 
4. 1 . Statement  of  the  Problem. 

A homogeneous,  laterally  unbounded  medium  is  assumed  to  occupy  the  region  z > 0.  A 
narrow  monoenergetic  pencil  beam  of  protons  is  assumed  to  be  incident  along  the  z-axis.  The 
information  sought  includes: 

(1)  dE/dz,  the  average  rate  of  energy  loss  by  the  proton  beam  per  unit  depth,  as  a 

function  of  the  depth  z.  The  units  of  dE/dz  are  MeV  cm^/g.  Actually  two 
such  quantities  are  calculated: 

(dE/dz),,,  pertaining  to  energy  losses  that  result  from  Coulomb  interactions  with 
atoms, 

(dE/dz)n,  pertaining  to  energy  losses  that  result  from  non-elastic  interactions  with 
nuclei. 

(2)  y(T,z),  the  proton  energy  spectrum,  expressed  in  terms  of  the  average  track  length 
per  unit  depth  and  per  unit  energy,  as  a function  of  the  proton  energy  T and  the 
depth  z.  The  units  are  MeV”^ 

(3)  f(p,z),  the  radial  distribution  of  the  energy  lost  by  the  proton  beam  (in  Coulomb 
interactions),  as  a function  of  the  distance  p from  z axis. 

Because  of  the  short  range  of  the  secondary  electrons,  it  is  a very  good  t^proximation  to 
consider  the  quantity  (dE/dz)^  equal  to  the  energy  deposition  per  unit  depth  due  to  Coulomb 
interactions,  (dD/dz)^.  As  will  be  shown  in  a later  report,  the  energy  deposition  per  unit  depth 
due  to  nuclear  interactions,  (dD/dz)„,  is  approximately  equal  to  0.6  (dE/dz)^.  The  discount 
factor  0.6  arises  mainly  from  the  fact  that  a large  fraction  of  the  energy  used  for  nuclear 
reactions  eventually  escapes  in  the  form  of  penetrating  secondary  neutrons. 

Also  computed  by  PTRAN3D  are  the  mean  energy  Tgv(2)  the  standard  deviation 
aj(z),  which  are  evaluated  as  track  averages: 


(4.1.1) 


(4.1.2) 
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4.2.  Track  Simulation. 


The  simulation  of  proton  tracks  is  based  on  the  procedures  outlined  in  Section  2,  and  uses 
the  data  sets  generated  by  programs  VPREP4  and  MPREP4.  Suppose  the  proton  has  reached 
energy  E,  and  the  events  in  the  n*^  step  are  to  be  sampled.  The  step  size  As^  is  obtained  by 
interpolating  to  energy  E in  the  table  of  fis^^  vs  T^j.  The  grid  index  m is  calculated  such  that  T^, 
is  the  closest  grid  energy  to  E.  The  scaled  energy  loss  X is  sampled  from  the  Vavilov 
distribution,  and  the  scaled  deflection  angle  from  the  Moli^re  distribution,  using  the  alias- 
method  arrays  for  grid  interval  (T^,  T^+j).  The  use  of  pre-computed  distributions  for  this  grid 
interval,  instead  of  a distribution  for  an  initial  energy  E and  path  length  As„,  is  an  ^proximation 
which  speeds  up  the  calculation.  The  error  incurred  thereby  is  very  small  because  the  variation 
with  energy  of  the  scaled  distributions  is  very  slow.  From  the  scaled  energy  loss  X,  the  actual 
energy  loss  E^  — En+j  is  determined  according  to  Eqs  (A2.1)  and  (A2.4)  of  Appendbc  2,  using 
A^v  = ” ^m+i»  values  of  ^ and  interpolated  to  energy  E.  The  multiple  scattering 

deflection  angle  6’  is  determined  from  i?  according  to  Eq  (A1.2)  of  Appendix  1,  using  a value  of 
Xc  interpolated  to  energy  E. 

4.3.  Crossing  of  Scoring  Planes. 

A set  of  scoring  planes  at  depths  z = z^,  f = 1,2,....,  is  introduced.  It  is  convenient  to 
define  these  depdis  in  terms  of  fractions  b^  of  the  CSDA  range  rQ  (at  the  initial  energy  of  the 
protons).  The  scaling  of  depths  in  units  of  rQ  minimizes  the  explicit  dependence  of  dE/dz  and 
y(T,z)  on  the  energy  of  the  incident  beam,  and  makes  it  easy  to  choose  scoring  planes  which 
allow  a good  description  of  the  Bragg  peak  of  the  depth-dose  distribution. 

The  values  of  b^  should  extend  from  zero  to  at  most  1.04.  Penetration  beyond  this  depth 
is  extremely  unlikely  and  of  no  practical  interest.  It  desirable  to  space  the  scoring  depths 
especially  closely  in  the  neighborhood  of  the  Bragg  peak,  which  occurs  for  b^  close  to  0.99. 

Suppose  that  the  n^  step  of  the  condensed  random  walk  involves  a crossing  of  the 
boundary  plane  z = Zf.  In  other  words,  z^  < Z;  < z^+i.  Let  E^  and  Ej,+i  denote  the 
corresponding  energies,  and  cos0„  and  cos0n^.2  Ae  corresponding  obliquity  cosines,  before  and 
after  the  crossing.  Let  a = (z^  - Zn)/(^n+i  ' ^n)  denote  the  fraction  of  the  step  that  has  been 
traversed  by  the  proton  when  the  crossing  occurs.  The  energy  E^^  at  the  point  of  crossing  is 
obtained  by  linear  interpolation, 

= E„  - a(E„-E„.,)  . (4.3.1) 


The  survival  factor  W^j.  = W(Egr),  the  nuclear  attenuation  factor  = fi(E^r),  and  the  stopping 
power  Sgj  = S(Egj)  are  obtained  by  interpolation  to  energy  E^^.  The  obliquity  cosine  at  the  point 
of  crossing  is  obtained  by  liner  interpolation, 

cos^gr  = cos^n  + a(cos0n+2  - cos^n)  . (4.3.2) 


The  radial  distance  at  the  point  of  crossing  is  calculated  as 


Per 


+ 


(Yn  ^ «(yn^l 


(4.3.3) 
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4,4.  Crossing  Scores. 

Four  scores  are  recorded  for  the  crossing,  which  depend  on  the  cosine  of  the  direction  of 
motion,  cos^g^,  and  on  the  nuclear  survival  weight  factor,  W^.^,  and  the  stopping  power 
evaluated  at  the  crossing  energy 


(1) 

Scr  Wer/cos0„ 

used  to  estimate  (dE/dz)^,  and  f(p,z); 

(2) 

^cr  PcT 

used  to  estimate  (dE/dz)^; 

(4.4.1) 

(3) 

1/cOS^cr 

used  to  estimate  y(T,z); 

(4) 

SJcosd^^ 

used  to  estimate  S(T)  y(T,z),  where  S(T)  is  the  stopping 
power. 

Actually  the  absolute  value  of  cos^^^  should  be  used  in  these  scores,  but  negative  values  of 
cos^gj.  were  never  encountered  in  all  the  runs  of  PTRAN  for  beam  energies  from  250  MeV  to 
50  MeV.  Radial  distributions  are  scored  by  assigning  Score  (1)  to  the  appropriate  radial 
histogram  bin  according  the  value  of  Energy  spectra  are  scored  by  assigning  Scores  (3)  and 
(4)  to  the  appropriate  energy  histogram  bin  according  to  the  value  of  Score  (4)  is  not  really 
necessary,  because  S(T)  y(T,z)  could  also  be  obtained  by  multiplying  the  Monte  Carlo  estimate  of 
y(T,z)  with  S(T)  at  the  end  of  the  calculation.  However,  the  use  of  (4)  is  inherently  more 
accurate. 

4.5.  Scoring  of  Track  Ends. 

When  the  proton  energy  E falls  below  a chosen  value  E^^j,  the  history  is  terminated,  and 
the  score  Eg^  W(Eg„t)  is  recorded  at  the  appropriate  depth  interval  of  a histogram.  It  has  been 
found  that  diis  histogram,  as  a function  of  depth,  can  be  represented  with  sufficient  accuracy  by  a 
Gaussian  distribution.  The  track-end  scores  are  utilized  by  a processing  program  called  PTSUM, 
which  is  discussed  in  Section  4.9. 


4.6.  Input  for  PTRAN3D. 


Input  Parameters.  The  user  is  prompted  to  supply  the  following  input  data: 


NBEG 

NFIN 

IHIST 


= Grid  index  that  specifies  the  initial  proton  energy.  The  relation  between  grid 
indices  and  energies  is  given  in  file  PARAM4.TBj  for  data  set  j, 

= Grid  index  for  the  fmal  proton  energy  at  which  the  Monte  Carlo  histories  are 
terminated 

= Number  of  Monte  Carlo  histories  to  be  sampled 


IMONIT  = Number  of  histories  in  group;  the  set  of  sampled  histories  is  divided  into 
(IHIST/IMONIT)  groups.  After  the  completion  of  each  successive  group, 
the  number  of  completed  histories  is  shown  on  the  monitor  screen.  Results 
for  different  groups  are  used  to  determine  the  statistical  uncertainties  of 
various  quantities. 
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INRAN  = Random  number  seed  for  congruential  random  number  generator  (should  be 
odd  in  order  not  to  reduce  the  period) 


BDFIL 


Name  of  the  boundary-information  file.  This  file  must  be  prepared  in 
advance,  and  should  contain  the  following  information: 


LBMAX,  IBMAX 


Number  of  scoring  depths  and  radial  bins 


(BR(L),  L=l, LBMAX) 


Depths  of  the  scoring  planes,  in  units  of 
the  CSDA  range 


(MEMAX(L),L= 1 ,LBMAX) 


Number  of  energy  bins  (of  equal  size) 
used  to  classify  proton  energy  spectra 


(EBOT(L),  L=l, LBMAX) 


Lowest  energies  (MeV)  used  for  spectral 
energy  classification 


(ETOP(L),  L=  1 , LBMAX) 
(RTOP(L),  L=l, LBMAX) 


Highest  energies  (MeV)  used  for  spectral 
energy  classification. 

Largest  radial  distance  (g/cm^)  for 
classifying  the  radial  distribution  of  energy 
loss. 


Cases  in  which  the  crossing  energy  is  smaller  than  EBOT  or  larger  than  ETOP  are  not 
included  in  the  spectrum.  Cases  in  which  the  radial  distance  at  the  time  of  crossing  exceeds 
RTOP  are  not  included  in  the  radial  distribution.  Appropriate  values  of  the  arrays  EBOT(L), 
ETOP(L)  and  RTOP(L)  can  be  estimated  from  the  results  of  preliminary  results  of  PTRAN3D. 
This  will  be  discussed  below  in  Section  4.8. 


Options.  The  user  is  prompted  by  the  program  to  decide  whether  energy-loss  straggling, 
angular  multiple-scattering  deflections,  and  spatial  multiple-scattering  displacements  should  be 
included.  Except  for  exploratory  model  studies,  the  answer  to  all  three  questions  should  be 
affirmative.  Finally,  the  user  is  asked  to  specify  the  name  of  the  output  file. 

Input  Arrays.  The  user  must  supply  the  names  of  three  large  input  files.  File  names  are 
suggested  by  the  program,  but  the  use  of  these  names  is  not  mandatory: 

PD  AT:  Contains  cross  sections  and  parameters,  generated  by  program  PARAM4; 

VDAT:  Contains  arrays  needed  for  sampling  from  the  Vavilov  energy-loss  distribution, 
generated  by  program  VPREP4; 

MDAT:  Contains  arrays  needed  for  sampling  from  the  Moli&re  multiple  scattering 
distribution,  generated  by  program  MPREP4. 

Also  available  must  be  two  auxiliary  input  files:  THSET2  contains  a set  of  reduced  angles 
required  for  the  sampling  of  angular  multiple  scattering  deflections.  GAUSS.DAT  contains  data 
used  for  selecting  random  variates  from  a Gaussian  distribution,  which  is  required  for  the 
sampling  of  lateral  multiple  scattering  displacements  according  to  Eqs  (2.7)  and  (2.8). 
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4.7,  Output  from  PTRAN3D, 

File  Names.  The  output  file  starts  with  a listing  of  the  names  of  the  input  files  used  for 
PDAT,  VDAT,  MDAT  and  BDFIL.  Then  follows  the  name  of  the  output  file.  With  this 
information  the  user  has  a complete  record  of  the  input  data  associated  with  each  output  file. 

Irq)ut  Parameters  and  Options.  Next  follows  a listing  of  various  input  parameters  and 
options  used  in  the  run,  including  IHIST,  INRAN,  LBMAX,  IBMAX,  NBEG,  NFIN  (described 
in  Section  4.6).  In  addition,  three  other  parameters  are  listed  which  indicate  what  options  were 
chosen: 


ISTRAG 

IMULT 

IREF 


1:  Energy-loss  straggling  included 
2:  Energy-loss  straggling  omitted 

1:  Multiple-scattering  deflections  included 
2:  Multiple-scattering  deflections  omitted 

1:  Lateral  multiple-scattering  displacements  included 
2:  Lateral  multiple-scattering  displacements  omitted 


Other  Pertinent  Quantities.  Also  listed  are: 


EBEG 

EFIN 

RANGE 


Initial  energy,  MeV,  corresponding  to  index  NBEG 
Final  energy,  MeV,  corresponding  to  index  NFIN 
CSDA  range,  g/cm^,  at  energy  EBEG. 


Output  Table.  The  next  section  of  the  output  file  has  the  form  of  a table  containing  the 
following  quantities: 


L 

BR 

SCORC  = 

PERCC  = 

SCORN  = 


Index  of  scoring  plane 

Depth  of  scoring  plane,  in  units  of  CSDA  range; 

(dE/dz)g,  mean  energy  loss  per  unit  depth,  MeV  cm^/g,  due  to  Coulomb 
interactions  with  atoms,  associated  with  crossing  of  the  L^  scoring  plane 

Error  (percent  standard  deviation)  of  SCORC,  estimated  from  the  dispersion 
of  the  results  obtained  with  (IHIST/IMONIT)  groups  of  histories 

(dE/dz)n,  mean  energy  loss  per  unit  depth,  MeV  cm^/g,  due  to  nonelastic 
nuclear  interactions,  associated  witii  the  crossing  of  the  L*^  scoring  plane 


PERCN  = Error  (percent  standard  deviation)  of  SCORN,  estimated  from  the  dispersion 
of  the  results  obtained  with  (IHIST/IMONIT)  groups  of  histories 

CROSS  = Fraction  of  proton  tracks  in  the  Monte  Carlo  simulation  which  crossed  the  L^ 
scoring  plane.  Note  that  nonelastic  interactions  do  not  prevent  such  crossings, 
because  they  are  taken  into  account  by  means  of  a survival  weight  factor. 


SCORC,  SCORN  and  CROSS  are  normalized  to  one  incident  proton. 
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RMAX 


RTOP 


EXCESS 


= Maximum  radial  distance  from  beam  axis  (g/cm^)  with  which  proton 
tracks  crossed  the  L’th  boundary  in  the  entire  Monte  Carlo  run. 
RMAX  is  also  stored  in  an  auxiliary  output  file  (see  Section  4.8).  If 
RTOP  is  properly  chosen,  RMAX  is  close  to  and  smaller  than 
RTOP. 

= Upper  limit  (g/cm^)  of  the  radial  distance  from  the  beam  axis 
included  in  the  histogram  representing  the  radial  distribution  of 
energy  loss 

= Fraction  of  the  energy  loss  at  the  scoring  depth  that  involves 
radial  distances  greater  than  RTOP.  Pertains  to  events  that  are  not 
included  in  the  radial  histogram. 


EBOT  (ETOP) 


Lower  (upper)  energy  limit  (MeV)  of  the  spectral  histogram 


EMIN  (EMAX)  = Minimum  (maximum)  energy  (MeV)  with  which  proton  tracks 

crossed  the  L’th  boundary  in  the  entire  Monte  Carlo  run.  EMIN  and 
EMAX  are  also  stored  in  an  auxiliary  output  file  (see  Section  4.8). 

If  EBOT  is  properly  chosen,  EMIN  is  close  to  and  larger  than 
EBOT,  and  EMAX  is  close  to  and  smaller  than  ETOP. 


MEMAX 


Number  of  spectral  energy  bins 


MLOW  (MHIG)  = Number  of  sampled  tracks  in  which  a crossing  energy  occurred  that 
was  smaller  than  EBOT  (or  larger  than  ETOP)  and  was  not  included 
in  the  computation  of  the  energy  spectrum 


Output  Arrays.  The  part  of  the  output  file  consists  of  large  arrays  for  use  in  other 
computer  calculations  that  analyze  the  output  of  PTRAN3D. 


(i)  Energy  Losses  in  Track  Ends 


DFSUM:  Energy, in  MeV  (per  incident  proton)  lost  in  track  ends 

KRMAX:  Number  of  depth  intervals  (of  equal  size)  in  which  the 

track-end  energy  losses  are  scored 

ZRMIN,  ZRMAX:  Minimum  and  maximum  depths  (in  units  of  the  CSDA 

range)  between  which  the  track-end  energy  losses  are 
scored 


(DUMPF(K),  K = l, KRMAX):  Track-end  energy  losses,  MeV,  in  various  depth  intervals 

DFSUM  and  DUMPF  are  normalized  to  one  incident  proton,  and  provide  input  needed  by 
the  processing  program  PTSUM  described  below  in  Section  4.6. 
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(ii)  Track  Lengths  and  Moments  of  Energy  Spectra 

(TRACK(L),  L=1,LBMAX)  = average  track  length  per  unit  depth  at  various  depths 

(EMOMl(L),  L=1,LBMAX)  = average  proton  energy,  MeV,  at  various  depths, 

calculated  as  track  average,  from  Eq  (4,1.1). 

(ESIG(L),  L=1,LBMAX)  = standard  deviation  of  the  proton  energy,  MeV,  at  various 

depths,  from  Eq  (4.1.2). 

TRACK  and  EMOMl  are  normalized  to  one  incident  proton.  In  the  computation  of 
EMOMl  and  ESIG  all  events  are  included,  even  the  rare  events  with  crossing  energies  E^^ 

smaller  than  EBOT  or  larger  than  ETOP. 

(iii)  Energy  Spectra  and  Radial  Distributions 

For  each  of  the  LBMAX  scoring  depths,  the  following  information  is  given: 


First  record:  L,  BR,  MEMAX,  EBOT,ETOP,  IBMAX,  RTOP,  where 


L 

BR 

MEMAX 
EBOT  (ETOP) 
IBMAX 
RTOP 


— index  of  scoring  plane 

— depth  of  scoring  plane,  in  units  of  CSDA  range 
= number  of  spectral  energy  bins 

= lower  (upper)  limit  of  spectral  histogram,  MeV 
= number  or  radial  bins 
= upper  limit  of  radial  histogram,  g/cnP’ 


Subsequent  records: 


(Y(M),  M = l, MEMAX) 


= histogram  of  y(T,z),  track  length,  per  unit  depth  and 
energy,  in  individual  energy  bins.  Units  are  MeV~^. 


(YS(M),  M = l, MEMAX)  = histogram  of  S(T)y(T,z),  track  length  multiplied  by 

stopping  power,  per  unit  depth  and  energy.  Units  are 
cnr/g. 


(RADC(I),I=1,IBMAX)  = fractions  of  the  energy  lost  in  each  of  the  IBMAX 

radial  bins.  Results  pertain  to  energy  losses  due  to 
Coulomb  interactions. 


The  histograms  Y(M)  and  YS(M)  are  normalized  to  one  incident  proton.  The  widths  of  the 
histogram  bins  is  DE  = (ETOP-EBOT)/MEMAX.  The  sum  Y(1)+Y(2)...+Y(MEMAX), 
multiplied  by  DE,  is  equal  the  track  length  per  unit  depth,  TRACK(L),  listed  previously  in  the 
output  file.  The  sum  YS(1)+YS(2)+...YS(MEMAX),  multiplied  by  DE,  is  equal  to  the 
contribution  to  (dE/dz)j.  from  scoring-plane  crossings. 
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The  histogram  RADC(I)  is  normalized  so  that  the  sum  RADC(1)+RADC(2)4-... 
+RADC(IBMAX)  + EXCESS(L)  is  equal  to  unity.  The  width  of  the  histogram  bins  is 
DR  = RTOP/IBMAX.  In  terms  of  the  radial  distribution  ftmction, 


(4.7.1) 


RADC(I) 


where  Rj  = (I-1)*DR  and  Rj  = I*DR. 

Timing.  The  last  two  items  of  the  output  file  consist  of: 

Elzqised  time  for  the  run,  in  minutes 

Dates  and  times  at  the  beginning  and  end  of  the  run. 

This  information  is  also  displayed  at  the  end  of  the  run  on  the  monitor  screen,  followed  by 
the  name  of  the  output  file. 

4.8.  Auxiliary  Output  File. 

An  auxiliary  output  file  is  automatically  written,  whose  name  is  the  same  as  that  of  the 
regular  output  file,  preceded  by  X.  This  file  contains  the  EMIN,  EMAX  and  RMAX  arrays  for 
all  the  crossing  boundaries.  TTie  arrays  obtained  in  a preliminary  run  of  PTRAN3D  can  be  used 
as  estimates  for  the  arrays  EBOT,ETOP  and  RTOP,  and  can  be  inserted  into  the  boundary 
crossing  file  BDFIL. 


4.9.  Program  PTSUM:  Combination  of  Crossing  and  Track  End  Scores 

A processing  program,  PTSUM,  uses  the  output  file  from  PTRAN3D  to  calculate  the  total 
energy  loss  (dE/dz)^  due  to  Coulomb  interactions  (practically  equal  to  energy  deposition  (dD/dz)^ 
by  combining  the  scores  from  scoring-boundary  crossings  and  from  track  ends.  PTSUM  also 
checks  the  accuracy  of  the  calculation,  by  verifying  that  (dE/dz)^  + (dE/dz)^,  integrated  over  all 
depths,  equals  the  energy  with  which  the  proton  enters  the  medium.  The  output  from  PTSUM  is 
stored  in  two  files:  a summary  table,  and  a file  with  arrays  useful  for  further  computer 
calculations.  Table  3 presents  a typical  summary  file,  for  the  case  of  a 160-MeV  proton  beam. 
The  contents  of  the  second  output  file  are  as  follows: 

Name  of  input  file  from  PTRAN3D 

Name  of  output  file  from  PTSUM 

Beam  energy  Tq  (MeV),  and  CSDA  range  rQ  (g/cm^)  at  energy  Tq 


UMAX 

(BR(L),L=1,LMAX) 
(SCORCT(L),L=l,LMAX)  = 
(SCORN(L),L=l,LMAX)  = 
(DF(L),L=1,LMAX) 


number  of  scoring  depths 

set  of  scoring  depths  (in  units  of  rg) 


set  of  values  of  (dE/dz)^ 
set  of  values  of  (dE/dz)^ 


contributions  to  (dE/dz)^  from  track  ends 
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Tables  4 and  5 list  results  obtained  with  PTSUM  for  proton  beams  with  seven  energies 
between  250  MeV  and  50  MeV.  Table  4 shows  the  fraaional  contribution  from  track  ends  to 
(dE/dz)g,  as  a function  of  depth,  and  Table  5 gives  the  ratio  (dE/dz)^  /(dE/dz)^  of  nuclear  to 
Coulomb  energy  losses,  again  as  a ftmction  of  depth. 


5.  Monte  Carlo  Program  PTRANID 

PTRANID  is  a one  dimensional  version  of  PTRAN3D,  which  omits  the  treatment  of  the 
radial  distribution  of  deposited  energy.  Again,  a narrow  pencil  beam  is  assumed  to  be  incident 
along  the  z-axis.  However,  because  of  the  omission  of  the  radial  variable  p,  results  obtained 
with  PTRANID  can  be  interpreted  as  applying  to  a broad  parallel  incident  beam. 

In  the  one-dimensional  treatment  involving  only  the  variable  z,  the  azimuthal  angle  (p  is 
irrelevant,  and  the  changes  of  azimuthal  angle  according  to  Eqs  (2.4)  and  (2.5)  need  not  be 
evaluated.  Moreover,  the  lateral  spatial  multiple-scattering  displacements  are  quite  unimportant. 
Therefore,  Eqs  (2.5)  and  (2.6)  are  not  used,  the  sampling  of  and  1^  from  a Gaussian 
distribution  is  avoided,  and  it  is  assumed  that  Ax'  — Ay'  ~ 0.  It  has  oeen  verified  that  the  error 
incurred  thereby  is  negligible  in  regard  to  (dE/dz)g,  (dE/dz)^  and  y(T,z).  As  a result  of  these 
simplifications,  runs  of  PTRANID  take  only  about  2/3  of  the  time  needed  with  PTRAN3D. 

PTRANID  requires  the  same  input  files  PD  AT,  VDAT  md  MDAT  as  PTRAN3D.  The 
boundary  input  file  BDFIL  can  also  be  the  same  as  that  used  with  PTRAN3D,  but  die  parameter 
IBMAX  and  the  array  RTOP,  pertaining  to  the  radial  distribution,  can  be  omitted.  TTie  auxiliary 
input  file  GAUSS.DAT  is  not  required. 

The  output  from  PTRANID  has  the  same  format  and  content  as  that  from  PTRAN3D, 
except  that  die  radial  energy-loss  distributions  RADC  are  omitted.  The  processing  program 
PTSUM  can  also  be  applied  to  the  output  from  PTRANID, 


6.  Random  Numbers 

Congruential  Generator.  The  random  numbers  needed  in  PTRAN3D  and  PTRANID  are 
obtained  by  a congruential  generator,  which  is  coded  in-line  and  requires  only  two  statements, 
and  is  therefore  very  fast: 

m =IAND(MASK,IR*MULT) 

R=RNORM*IR 

IR  is  a random  integer.  R is  the  corresponding  real  random  number  (between  0 and  1) 
from  a uniform  distribution,  obtained  with  normalizing  factor  RNORM  = 2’^^ 

= 4.656612873E-10.  MASK  has  the  decimal  value  2147483647,  which  corresponds  to  a binary 
number  0111  1111  1111  1111  1111  1111  1111  1111.  The  operator  LAND,  ^plied  to  the 
product  IR*MULT,  makes  the  most  significant  bit  of  IR  zero.  The  multiplier  MULT  should  be 
chosen  so  that  there  is  a minimum  amount  of  correlation  between  successive  numbers  in  the 
pseudo-random  number  sequence.  The  period  of  the  sequence  can  be  shown  to  be 
2^^  = 5.37E-f-08,  The  initial  value  of  IR  (the  random  number  seed)  should  be  an  odd  integer; 
otherwise,  the  period  will  be  shortened. 
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Tests  with  program  PTRANID,  using  step-size  Grid  4,  indicated  that  approximately  1500 
random  numbers  were  required  for  sampling  a single  proton  track  from  an  initial  energy  of  160 
MeV  down  to  140  keV.  Thus  the  period  of  the  random  number  generator  is  such  that  only 
358,000  histories  can  be  generated  before  the  sequence  of  random  numbers  repeats  itself. 

In  PTRAN3D  and  PTRANID,  two  sequences  of  random  numbers  are  actually  used:  IRA 
(with  a multiplier  MULTA  = 69069)  and  IRB  (with  a multiplier  MULTB  = 1664525). 

According  to  Knuth  (1981)  these  multipliers  give  rise  to  sequences  which  have  a satisfactory  lack 
of  correlation,  as  indicated  by  the  theoretical  spectral  test.  An  IRB  sequence  is  used  to  provide 
random  number  seeds  for  successive  proton  histories.  Within  each  history  random  numbers 
from  an  IRA  sequence  are  used. 

With  this  procedure  one  can  generate  approximately  537  million  histories  each  of  which  has 
a different  starting  random  number  seed.  One  cannot  guarantee  that  there  are  no  repetitions  and 
overlaps  between  the  random-numbers  sets  for  different  histories  generated  with  the  multiplier 
MULTA.  It  seems  very  unlikely,  however,  that  such  repetitions  will  occur  often  enough  so  that 
they  could  significantly  influence  the  Monte  Carlo  results.  The  best  way  to  reassure  oneself  is  to 
m^e  comparisons  with  calculations  in  which  random  numbers  are  generated  by  a different 
generator  with  a longer  period. 

Fibonacci  Generator.  An  alternative  version  of  PTRANID,  called  PTRANIDX,  uses  a 
random  number  generator  due  to  Marsaglia  and  Zaman  (1987).  This  generator  is  based  on  a 
Fibonacci  sequence,  and  is  stated  to  satisfy  stringent  empirical  tests  of  randomness.  The 
generator  has  the  very  long  period  2^  = 1.79E13.  In  PTRANIDX  a subroutine  RANMAR  from 
James  (1988)  is  used,  which  is  an  adaptation  of  the  Marsaglia-Zaman  generator  and  generates 
random  numbers  in  batches  of  1000,  thereby  reducing  the  time  expended  on  calls  to  the 
subroutine.  There  is  also  an  initializing  subroutine  RMARIN  which  requires  two  input  numbers 
as  seeds. 

Several  numerical  experiments  were  carried  out  comparing  results  for  the  same  problem 
obtained  with  PTRANID  and  PTRANIDX.  No  statistically  significant  differences  were  found  in 
regard  to  (dE/dz)^,  (dE/dz)„,  and  y(T,z).  However,  the  required  computing  time  was  about  1/3 
higher  for  PTRANIDX  than  for  PTRANID. 

Use  of  Intrinsic  Random-Number  Generating  Function.  Many  Fortran  compilers  provide  an 
intrinsic  function  that  generates  uniform  random  numbers,  usually  by  the  multiplicative 
congruential  method.  As  a convenience  for  users  who  would  like  to  use  such  an  intrinsic 
function,  another  version  of  PTRANID,  called  PTRANIDY,  was  prepared.  In  this  version 
random  numbers  are  assumed  to  be  generated  by  an  intrinsic  function  RND,  and  the  random 
number  sequence  is  started  by  an  intrinsic  function  RANDS  (x),  where  x (a  real  number  between 
0 and  1)  is  the  random  number  seed. 
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7.  Influence  of  Step-Size  Grid 


Runs  of  program  PTRANID,  for  a proton  beam  energy  with  an  initial  energy  of  160  MeV, 
were  made  with  five  step-size  grids,  as  described  in  Section  3.1.  The  grid  parameter  AT  and  the 
computing  times  for  sampling  (and  scoring)  1 million  Monte  Carlo  histories  down  to  an  energy  of 
140  keV  was  as  follows: 


Grid 

AT 

(MeV) 

Time 

(Minutes) 

Relative 

Time 

1 

4 

162 

1.0 

2 

2 

190 

1.17 

3 

1 

263 

1.62 

4 

0.5 

475 

2.93 

5 

0.25 

750 

4.63 

These  runs  were  made  using  a personal  computer  with  a 25  Mhz  486  processor  and  a Weitek 
4167  coprocessor. 

In  regard  to  the  quantity  (dE/dz)^  the  differences  between  the  results  obtained  with  different 
grids  were  small:  less  than  0.1  percent  for  depths  z < 0.9  ro;  less  than  1 percent  for  0.9 
ro  < z < 1.01  rQ,  and  less  than  3 percent  for  z = 1.02  rQ. 

Larger  differences  were  found  in  regard  to  the  track-length  distribution  y(T,z).  Figs.  3,  4 
and  5 , for  depths  z = 0.1  rQ,  0.5  rQ,  and  0.9  rQ,  respectively,  for  a beam  energy  of  160  MeV, 
shows  the  percentage  amounts  by  which  y(T,z)  obtained  with  Grid  2,  Grid  3,  or  Grid  4 differs 
from  the  value  obtained  with  Grid  5.  The  differences  are  1 to  3 percent  in  the  central  region  of 
the  spectrum,  but  are  considerably  larger  in  the  wings  of  the  spectrum. 

Comparison  runs  were  also  made  of  radial  distributions  calculated  with  program 
PTRAN3D,  with  Grids  1,  2,  3,  and  4.  Results  are  shown  in  Figs. 6 and  7 for  a 160-MeV  beam 
of  protons,  at  a depth  z = 0.1  rQ.  In  addition  to  the  Monte  Carlo  histograms,  smooth  curves  are 
shown  which  represent  the  radi^  distribution  calculated  from  the  theory  of  Moli^re  (1955),  with 
energy-loss  taken  into  account  in  the  continuous-slowing-down  approximation. 

The  Monte  Carlo  results  in  figure  6 were  obtained  without  taking  into  account  the  spatial 
multiple-scattering  displacements,  i.e.,  setting  Ax'  = Ay'  = 0 in  Eq  (2.10).  It  can  be  seen  that 
the  finer  the  grid  the  closer  the  histograms  approach  the  curve  from  the  Moli&re  theory.  In 
figure  7 a similar  set  of  results  is  shown,  which  were  obtained  with  the  lateral  multiple-scattering 
displacements  included.  In  this  case  there  is  little  variation  from  one  grid  to  another. 

It  is  interesting  that  the  Monte  Carlo  results  obtained  with  Grid  4 (in  fig.  6)  or  with  all 
grids  (in  fig.  7),  agree  well  with  the  curve  from  the  Moli^re  theory.  Similar  agreement  was  also 
obtained  at  other  depths.  The  implication  is  that  the  use  of  PTRAN3D  is  not  really  necessary, 
but  could  be  replaced  by  the  combined  use  of  PTRANID  and  the  Moli^re  theory.  Of  course  this 
could  only  be  done  for  depths  smaller  than  the  CSDA  range.  Beyond  the  Bragg  peak  the 
PTRAN3D  program  is  still  needed. 
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8.  Program  and  Data  Files 


Users  can  obtain  from  the  author  of  this  report  a 3.5  inch,  1.44  Mb  disk  with  four  archive 
flies,  named  PARAM.EXE,  VPREP.EXE,  MPREP.EXE  and  PTRAN.EXE.  The  contents  of 
these  files  are  listed  in  table  7.  The  archive  files  are  self-extracting.  For  example,  by  running  the 
program  PARAM  one  recovers  all  the  individual  files  stored  in  the  archive  file  PARAM.EXE. 

Included  among  the  files  in  VPREP.EXE  and  MPEP.EXE  are  the  formatted  files 
VPREP4.ARj  and  MPREP4.ARj,  j = 2,3,4,  but  not  the  corresponding  binary  files  VPREP4.BRj 
and  MPREP4.BRj  which  are  needed  to  run  PTRANID  or  PTRAN3D.  The  reason  for  this 
omission  is  that  binary  files  are  in  general  not  portable,  and  depend  on  the  type  of  computer,  and 
on  the  compiler  used.  The  user  can  easily  generate  the  binary  files  with  the  VCON4  and 
MCON4  programs. 
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Appendix  1.  Moli^re  Distribution 


According  to  the  theory  of  Molibre  (1948),  the  distribution  of  multiple-scattering  angular 
deflections  can  be  expressed  as  a unique  function  of  a scaled  angular  variable  t?: 


FM(^)t5di? 


t?dt? 


... 

B 


(ALl) 


The  relation  between  the  scaled  angle  and  the  actual  multiple-scattering  deflection  angle  6 is 


(A1.2) 


where  Xc  and  B are  functions  of  the  particle  energy  and  the  path  lengdi  As. 

Moli^re’s  theoiy  was  developed  in  the  small-angle  approximation,  and  is  applicable  when 
the  multiple  scattering  angles  6 are  no  greater  than  approximately  20  degrees.  In  the  present 
application  to  protons  this  is  not  a significant  restriction.  The  expansion  in  inverse  powers  of  B 
in  Eq  (Al.l)  is  accurate  only  when  B is  larger  than  about  4.5.  As  discussed  by  Molifere,  the  use 
of  additional  terms  in  the  expansion  would  not  significantly  increase  the  accuracy,  because  of 
^proximations  made  elsewhere  in  the  derivation  of  theory. 

The  expansion  coefficients  in  Eq  (Al.l)  are  given  by 


i I “ ydy  Jo(i?y) 


(A1.3) 


when  Jq  denotes  a Bessel  function. 

The  quantities  f^^^  and  f^^^  were  first  tabulated  by  Moli^re  (1948),  and  later  by  Betiie  (1953) 
and  Scott  (1963).  In  the  course  of  the  present  work  an  even  more  extensive  table  of  coefficients 
was  calculated.  This  table  includes  values  of  f^^^  and  f^^^  at  131  values  of  ^ between  0 and  2.6, 
and  331  values  of  and  at  331  values  of  between  2.4  and  60.0.  The  new  results 
are  contained  in  the  files  MOLCl.COF  and  MOLC2.COF  distributed  with  the  PTRAN 
programs. 


For  a compound,  Xc  is  obtained  as  a sum  over  the  corresponding  quantities  for  the  atomic 
constituents: 

2 2 


Xc  ~ ^cj  ’ 


(A1.4) 
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where  Wj  is  the  fraction  by  weight  of  the  j**'  atomic  constituent,  and 

4 = 

Here  is  the  Avogadro  constant,  and  Aj  are  the  atomic  number  and  atomic  weight  of  the 
constituent,  r^  is  the  classical  electron  radius,  m/M  is  the  electron-proton  mass  ratio,  r is  the 
proton  kinetic  energy  in  units  of  the  proton  rest  mass,  and  As  is  the  path  length. 

The  Moli^re  parameter  B depends  on  the  ratio  of  the  characteristic  angle  Xc  ^ ^ screening 
angle  Xa>  is  obtained  from  the  equation 

B -logB  = log  ixJXaf’  + 1-  27  , (A1.6) 


m 

M 


T-t-1  1 
t(t+2)J 


a 

Aj 


As 


(A1.5) 


where  7 = 0.5772156649., , is  Euler’s  constant.  This  equation  can  be  easily  solved  iteratively  by 
Newton’s  method. 


The  screening  angle  Xa  is  given  by 


with 


logx^  = (l/x3  E WjXej  [logx^j  - Fj/Zj 


‘■aj 


■ 

2 

m 

a 

M 

kxF 

1.13+3.76  (ZjO//3)^ 


Zi 


2/3 


t(t  + 2) 


(A1.7) 


(A1.8) 


Here  k^p  = (971^)^'^^  2”^^^  = 0.88534...  is  a constant  associated  with  the  Thomas-Fermi  model, 
a is  the  fine-structure  constant  and  /3  is  the  proton  velocity  in  units  of  the  velocity  of  light.  The 
quantity  in  square  brackets  in  Eq  (A1.8)  is  an  approximation  obtained  by  Moli^re  with  a 
Thomas-Fermi  potential.  The  correction  factor  kjjp  (which  depends  on  Z^aiP)  converts  the 
screening  angle  to  one  corresponding  to  a Hartree-Fock  potential,  and  was  calculated  by  Berger 
and  Wang  (1988),  using  a formulation  of  Molifere’s  theory  given  by  Zeitler  and  Olsen  (1964). 
Values  of  kjjp  for  selected  elements  are  given  in  table  8. 

The  term  Fj/Zj  in  Eq  (A1.7)  is  an  addition  to  Moli6re’s  theory,  due  to  Fano  (1956),  that 
takes  into  influence  of  the  orbital  electrons.  Fj  is  given  by 


Fj  = log[ll30/?2  (1-^^)'^ 


(A1.6) 


where  the  constant  Uj  has  value  -3.6  for  hydrogen  and  -5.1  for  oxygen. 
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Appendix  2.  Vavilov  Distribution 


The  distribution  of  energy  losses  from  the  theory  of  Vavilov  (1957),  Fy(X,/3,/c),  is  a 
function  of  a scaled  energy-loss  variable  X,  and  also  depends  — more  weakly  — on  0 and  on  a 
skewness  parameter  k defined  below.  For  small  k,  the  Vavilov  distribution  ^proaches  the 
Landau  distribution  (Landau,  1944),  whereas  for  large  k it  approaches  a Gaussian  distribution. 

The  relation  between  the  scaled  energy-loss  X and  die  energy  loss  A is 


X = 


A-A 


av 


+ X, 


av 


(A2.1) 


where  A^y  is  the  average  energy  loss,  and 

{ , (A2.2) 

j 

Here  r^  is  the  classical  electron  radius,  mc^  the  electron  rest  energy,  and  Aj,  and  Wj  are  the 
atomic  number,  atomic  weight  and  fraction  by  weight  of  the  j*  atomic  constituent.  The  quantity 

Xgy  in  Eq  (A2.1)  is 


^av  = + 7 - 0^ 


- log  /c  , 


(A2.3) 


where  y = 0.5772156649...  is  Euler’s  constant.  The  skewness  par^neter  k is 


where 


-M 


2mc^/3^ 

1-/3^ 


9 


(A2A) 


(A2.5) 


is  the  maximum  amount  of  energy  which  a proton  can  lose  in  a collision  with  an  orbital  electron 
(considered  as  free). 

The  quantity  ^ has  the  following  significance.  In  a path  length  As,  diere  will  occur,  on  the 
average,  one  collision  with  an  orbital  electron  in  which  the  proton  loses  an  amount  of  energy 
greater  than  The  Vavilov  distribution  is  accurate  under  two  conditions:  (1)  the  path  length 
must  be  long  enough  so  that  | is  much  greater  than  the  mean  excitation  energy  of  the  material 
(75  eV  for  water);  (2)  The  path  length  must  be  short  enough  so  that  the  me^  energy  loss  A^y  is 
small  compared  to  the  initial  proton  energy. 
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The  formula  for  the  Vavilov  distribution  is 


Fv(X,/3,K)dX  = -^  j expjhj  (y)]  cos[h2  (y)]  dy  , (A2.6) 

where 


hiO-)'*  fi(y)  - 

» 

(A2.7) 

h2(y)  = «[y(X  + log/c)  + f2(y)]  , 

(A2.8) 

fi(y)  = /3^[logy  - Ci(y)]  - cosy 

-ySi(y)  , 

(A2.9) 

= y [log  y “ Ci(y)]  + sin  y + 

^^Si(y)  , 

(A2.10) 

Si(y)  = f ^ dt  (sine  integral)  , 

Jo  t 

(A2.11) 

f ^ cos  t 

Ci(y)  = 7 + logy  - dt  (cosine  integral)  . 

Jo  t 

(A2.12) 

Efficient  methods  for  evaluating  the  sine  and  cosine  integrals  by  series  expansions  and  rational 
^proximations  can  be  found  in  Abramowitz  and  Stegun  (1964). 

The  quantity  € in  Eq  (A2.7)  is  a binding  correction  which  was  applied  Blunck  and 
Leisegang  (1951)  to  Landau’s  theory  and  by  Shulek  et  al.  (1959)  to  Vavilov’s  theory.  As  shown 
by  Fano  (1963)  the  variance,  0^,  of  the  energy-loss  distribution  is 

={Em(1-^^/2)(1*«)  , (A2-‘3) 


with 


2S, 


e = 


Em(1  -r/2) 


log 


2mc 


(A2.14) 
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Sj  and  Ij  are  averages  taken  over  the  oscillator  strength  distribution: 


(A2.15) 

(A2.16) 


where  df/dE  is  the  density  of  optical  dipole  oscillator  strength  per  unit  excitation  energy,  E, 
above  the  ground  state.  In  the  present  work,  experimental  values  for  these  quantities  were  used, 
from  a compilation  of  Zeiss  et  al.  (1977):  Sj  = 191.2  eV,  Ij  = 939.9  eV.  These  values  were 
obtained  for  water  vapor,  but  is  expected  that  their  use  for  liquid  water  will  cause  little  error. 
With  these  input  data  one  obtains  the  following  values  of  e: 

Energy  (MeV)  = 200  100  50  20  10  5.13  3.97 

€ = 0.00554  0.00973  0.0168  0.0238  0.0553  0.0849  0.0984 

Straggling  measurements  by  Besenbacher  et  al.  (1981)  for  protons  in  low-Z  gases  suggest 
that  at  low  energies  Eq  (A2.14)  overestimates  12^.  In  the  present  work,  Eq  (A2.14)  was  therefore 
used  only  down  to  an  energy  where  e reaches  a limiting  value  = 0.1.  Below  that  energy,  e 
was  kept  constant. 
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Table  1.  Proton  stopping  powers  and  ranges  in  water.  Calculated  assuming  a mean  excitation 
energy  I = 75  eV. 


T - Proton  energy,  MeV 

STOP(e)  - electronic  stopping  power,  MeV  cmVg 
STOP(n)  - nuclear  stopping  power,  MeV  cmVg 
STOP(t)  - total  stopping  power,  MeV  cm^/g 
RANGE (c)  - CSDA  range,  cmVg 
RANGE (p)  - projected  range,  cm^/g 

DETOUR  - detour  factor  - RANGE(c)/RANGE(p) 


T 

STOP(e) 

STOP(n) 

STOP(t) 

RANGE(c) 

RANGE (p) 

DETOUR 

0.1000 

8 . 145E+02 

1.620E+00 

8.161E+02 

1.607E-04 

1.458E-04 

0.9073 

0.2000 

6.604E+02 

9.016E-01 

6.613E+02 

2.966E-04 

2.806E-04 

0.9460 

0.3000 

5.497E+02 

6.351E-01 

5 . 504E+02 

4.631E-04 

4.462E-04 

0.9635 

0.4000 

4.714E+02 

4.928E-01 

4.719E+02 

6.599E-04 

6.422E-04 

0.9731 

0.5000 

4.128E+02 

4.043E-01 

4.132E+02 

8.869E-04 

8.683E-04 

0.9790 

0.6000 

3.676E+02 

3.438E-01 

3.680E+02 

1.144E-03 

1.124E-03 

0.9829 

0.8000 

3.037E+02 

2.658E-01 

3.039E+02 

1.745E-03 

1.724E-03 

0.9877 

1.0000 

2.606E+02 

2.173E-01 

2 . 608E+02 

2.458E-03 

2.435E-03 

0.9905 

2.0000 

1.585E+02 

1.157E-01 

1.586E+02 

7.555E-03 

7.519E-03 

0.9952 

3 . 0000 

1.171E+02 

7.972E-02 

1.172E+02 

1.499E-02 

1.494E-02 

0.9965 

4.0000 

9.398E+01 

6.113E-02 

9.404E+01 

2.458E-02 

2.451E-02 

0.9971 

5 . 0000 

7.906E+01 

4.970E-02 

7.911E+01 

3.623E-02 

3.613E-02 

0.9974 

6.0000 

6.854E+01 

4.195E-02 

6.858E+01 

4.984E-02 

4.972E-02 

0.9976 

8.0000 

5.456E+01 

3.208E-02 

5.460E+01 

8.277E-02 

8.259E-02 

0.9978 

10.0000 

4.564E+01 

2.603E-02 

4.567E+01 

1.230E-01 

1.228E-01 

0.9980 

15.0000 

3.290E+01 

1.778E-02 

3.292E+01 

2.539E-01 

2.535E-01 

0.9982 

20.0000 

2.605E+01 

1.356E-02 

2.607E+01 

4.260E-01 

4.252E-01 

0.9983 

30.0000 

1.875E+01 

9.239E-03 

1.876E+01 

8.853E-01 

8.839E-01 

0.9984 

40.0000 

1.487E+01 

7.034E-03 

1.488E+01 

1.489E+00 

1.486E+00 

0.9985 

50.0000 

1.244E+01 

5.691E-03 

1.245E+01 

2.227E+00 

2.224E+00 

0.9985 

60.0000 

1.078E+01 

4.786E-03 

1.078E+01 

3.093E+00 

3.089E+00 

0.9986 

70.0000 

9.555E+00 

4.134E-03 

9.559E+00 

4.080E+00 

4.075E+00 

0.9986 

80.0000 

8.622E+00 

3.641E-03 

8.625E+00 

5 . 184E+00 

5.176E+00 

0.9986 

90.0000 

7.884E+00 

3.255E-03 

7.888E+00 

6.398E+00 

6.389E+00 

0.9986 

100.0000 

7.286E+00 

2.944E-03 

7.289E+00 

7.718E+00 

7.707E+00 

0.9987 

110.0000 

6.791E+00 

2.689E-03 

6.794E+00 

9 . 140E+00 

9.128E+00 

0.9987 

120.0000 

6 . 374E+00 

2.475E-03 

6.377E+00 

1.066E+01 

1.065E+01 

0.9987 

130.0000 

6.018E+00 

2.294E-03 

6.021E+00 

1.228E+01 

1.226E+01 

0.9987 

140.0000 

5.711E+00 

2.137E-03 

5.713E+00 

1.398E+01 

1.396E+01 

0.9987 

150.0000 

5.443E+00 

2.001E-03 

5.445E+00 

1.577E+01 

1.576E+01 

0.9987 

160.0000 

5 . 207E+00 

1.882E-03 

5.209E+00 

1.765E+01 

1.763E+01 

0.9988 

170.0000 

4.997E+00 

1.777E-03 

4.999E+00 

1.961E+01 

1.959E+01 

0.9988 

180.0000 

4.810E+00 

1.683E-03 

4.812E+00 

2.165E+01 

2.163E+01 

0.9988 

190.0000 

4.642E+00 

1.598E-03 

4 . 644E+00 

2.377E+01 

2.374E+01 

0.9988 

200.0000 

4.491E+00 

1.522E-03 

4.492E+00 

2.596E+01 

2.593E+01 

0.9988 

210.0000 

4.353E+00 

1.453E-03 

4.354E+00 

2.822E+01 

2.819E+01 

0.9988 

220.0000 

4.227E+00 

1.390E-03 

4.229E+00 

3.055E+01 

3.052E+01 

0.9988 

230.0000 

4.112E+00 

1.332E-03 

4.114E+00 

3.295E+01 

3.291E+01 

0.9988 

240.0000 

4.007E+00 

1.279E-03 

4.008E+00 

3 . 541E+01 

3.537E+01 

0.9988 

250.0000 

3.910E+00 

1.231E-03 

3.911E+00 

3.794E+01 

3.790E+01 

0.9989 

27 


V 

o 

X 

U 


C3 

E- 


iH 

rH 

iH 

rH 

fH 

fH 

fH 

iH 

fH 

fH 

O 

O 

O 

O 

o 

O 

O 

O 

O 

O 

O 

O 

> 

w 

M 

w 

w 

M 

w 

w 

w 

w 

w 

w 

w 

2 

O) 

O 

rs. 

iH 

CO 

o 

CO 

rs. 

fs. 

rs. 

rs. 

^S 

5 

40 

O 

o 

a> 

CD 

to 

to 

CD 

CO 

CO 

rs. 

CO 

00 

eg 

O) 

O) 

40 

eg 

eg 

N 

eg 

eg 

fs. 

0) 

40 

eg 

40 

CD 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

rs. 

o> 

CO 

00 

CO 

CD 

CD 

CD 

to 

CD 

CD 

o> 

00 

rs 

fSs 

CD 

CO 

CD 

CD 

to 

to 

CD 

CD 

CO 

CO 

eg 

eg 

eg 

eg 

eg 

o 

o 

o 

o 

o 

eg 

eg 

eg 

rg 

rg 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

+ 

+ 

+ 

+ 

+■ 

+ 

+ 

u 

w 

M 

M 

M 

M 

w 

w 

w 

w 

M 

M 

w 

< 

M 

w 

Ui 

Ui 

Ui 

Ui 

Ui 

Ui 

Ui 

IP 

IP 

IP 

CO 

rs. 

(O 

o 

CD 

CO 

CO 

o 

o 

o 

o 

o 

pg 

OO 

cn 

Oi 

CD 

00 

eg 

A 

4 

A 

o> 

o 

CO 

rv 

CD 

eg 

rss 

4 

CO 

o 

o 

o 

o 

o 

rs. 

cn 

eo 

fs. 

eg 

rg 

A 

cn 

eo 

O) 

4 

eg 

H 

O) 

rs. 

eg 

4 

4 

CO 

00 

o 

o 

o 

o 

o 

00 

eo 

eo 

QO 

Oi 

o 

Oi 

r^ 

CO 

CO 

4 

eg 

H 

f-* 

CO 

o 

OO 

rs» 

CD 

4 

o 

o 

o 

o 

o 

CD 

CD 

CO 

o 

a> 

Oi 

eg 

CO 

4 

eo 

a> 

< 

m 

rsk 

<3 

o 

CO 

r^ 

eg 

o 

o 

o 

o 

o 

4 

CD 

00 

CO 

00 

A 

eg 

eg 

eo 

4 

CO 

O) 

O) 

fH 

fH 

rg 

o 

o 

o 

o 

o 

<«• 

CO 

CO 

(O 

CO 

CO 

CO 

CO 

CO 

(O 

eo 

iH 

fH 

rg 

fH 

rg 

eg 

eg 

eg 

eo 

eo 

CO 

4 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

O 

o 

O 

o 

o 

o 

o 

o 

o 

o 

O 

w 

M 

M 

u 

M 

M 

M 

M 

w 

M 

w 

M 

M 

Ui 

Ui 

Ui 

w 

Ui 

M 

w 

IP 

IP 

IP 

IP 

W 

eg 

<4* 

OO 

40 

O) 

40 

eg 

N 

o 

to 

CO 

CO 

fs. 

eo 

Oi 

CD 

eo 

4 

rg 

4 

Oi 

4 

eo 

fs 

eg 

C3 

rs. 

o 

CD 

CO 

oo 

CO 

4 

00 

00 

CO 

iH 

fH 

QO 

4 

00 

<o 

o 

CD 

fs. 

eg 

fs. 

< 

«r 

40 

00 

4 

O 

CO 

00 

O) 

rv 

eg 

o 

CO 

r*^ 

eg 

eg 

Oi 

fH 

o 

Oi 

eg 

4 

CD 

rH 

fH 

C) 

oo 

iH 

CD 

40 

eg 

eo 

CO 

eg 

oo 

to 

CD 

rg 

00 

cn 

N 

rg 

4 

A 

fH 

4 

tH 

fH 

4 

CO 

CO 

O 

O 

CO 

o 

4 

fH 

eg 

eg 

eg 

iH 

0) 

CO 

4 

eo 

eg 

iH 

4 

eg 

4 

eg 

4 

r*. 

o> 

fH 

eg 

CO 

CO 

40 

CO 

lO 

CO 

4 

rH 

fH 

fH 

fH 

fH 

fH 

eg 

N 

N 

eg 

eg 

eg 

eg 

eg 

eg 

eg 

o 

O 

o 

O 

o 

o 

o 

O 

o 

O 

O 

O 

o 

O 

O 

o 

o 

o 

o 

o 

o 

o 

O 

O 

1 

H 

+ 

+ 

+ 

+ 

g- 

+ 

w 

u 

M 

M 

w 

w 

w 

w 

M 

u 

U4 

< 

Ui 

M 

Ui 

u 

Ui 

Ui 

Ui 

Ui 

IP 

IP 

IP 

IP 

40 

O) 

4* 

eg 

fH 

o 

O) 

N 

4 

4 

eg 

4 

2 

CO 

rg 

4 

eg 

cn 

o 

4 

o 

CO 

00 

o 

4 

< 

o 

40 

CO 

o 

CO 

eg 

fO 

fH 

to 

4 

4 

IH 

r^ 

4 

eg 

cn 

CD 

4 

eg 

G) 

eg 

4 

CD 

u 

o 

CO 

CO 

o 

fH 

eg 

CO 

CD 

to 

eg 

OO 

CO 

CO 

cn 

4 

CD 

4 

N 

rg 

o 

fs, 

o 

4 

fH 

CO 

CO 

CO 

fH 

O 

o 

r^ 

eg 

rg 

fH 

CO 

4 

X 

CO 

cn 

o 

eg 

CD 

eg 

00 

4 

to 

o 

4 

O 

CO 

40 

O 

o 

fH 

CO 

CO 

o 

rg 

rH 

CO 

CO 

eo 

4 

4 

4 

cn 

cn 

CO 

fH 

iH 

CD 

4 

rN 

eg 

CO 

40 

rs. 

rs. 

CD 

QO 

OO 

CD 

iH 

fH 

fH 

fH 

fH 

fH 

o 

o 

o 

O 

o 

eg 

eg 

eg 

eg 

N 

eg 

eg 

eg 

eg 

CO 

eo 

eo 

o 

O 

O 

O 

o 

O 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

1 

1 

w 

M 

w 

w 

w 

U4 

u 

w 

M 

w 

M 

Ui 

IH 

Ui 

Ui 

M 

Ui 

M 

Ui 

Ui 

Ui 

IP 

IP 

IP 

IP 

iH 

CO 

o 

CO 

CD 

O) 

CO 

40 

o 

o 

4 

Oi 

CO 

eg 

00 

Oi 

CD 

in 

CO 

in 

Oi 

fs. 

o 

00 

n 

eg 

CD 

CD 

CD 

rs» 

rs. 

00 

eg 

40 

eg 

o 

X 

eo 

CD 

CD 

CO 

00 

o> 

o 

Oi 

o 

eo 

CO 

•«T 

eg 

o 

rs. 

O) 

CO 

CD 

00 

4 

CO 

00 

Oi 

cn 

4 

oo 

cn 

fH 

CD 

cn 

m 

eg 

eo 

rH 

CO 

CO 

eo 

eg 

eg 

o 

O) 

40 

QO 

rg 

CO 

rs. 

oo 

Oi 

O 

fH 

4 

Oi 

eo 

m 

eg 

cn 

CO 

o 

fH 

fH 

o 

0) 

CO 

o 

CD 

rs. 

eg 

eg 

eo 

eo 

CO 

eo 

4 

eg 

fH 

fs. 

4 

eo 

iH 

fH 

fH 

fH 

fH 

fH 

O) 

CD 

rs. 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

eg 

eg 

fH 

fH 

CO 

eo 

eo 

CO 

eo 

eg 

eg 

eg 

eg 

eg 

eg 

eg 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

4* 

+ 

1 

1 

w 

M 

w 

Ui 

M 

44 

u 

w 

w 

W 

U 

u 

eg 

Ui 

w 

Ui 

Ui 

Ui 

M 

Ui 

Ui 

IP 

IP 

IP 

IP 

o 

o 

Oi 

C3> 

4 

o 

0) 

CM 

CD 

o 

fH 

eg 

eg 

H 

eg 

o 

fH 

Oi 

o 

O 

eg 

Oi 

eg 

fs. 

eo 

G> 

u 

40 

40 

fH 

CO 

rx 

4 

G) 

eg 

CO 

fH 

00 

A 

to 

cn 

CD 

O) 

O) 

to 

00 

eo 

o> 

4 

fs. 

CO 

CO 

Oi 

eg 

CD 

eg 

CO 

CO 

«H 

CO 

rH 

CD 

o 

W 

CO 

eg 

cn 

A 

o 

cn 

eo 

eg 

00 

Oi 

Oi 

Oi 

r*. 

4* 

o> 

(O 

4 

4 

O 

rs. 

fH 

CO 

4 

CO 

a 

oo 

eo 

o 

4 

o 

CD 

rs. 

eg 

00 

ca 

Oi 

Oi 

CD 

40 

CO 

rs. 

0) 

o 

4 

QO 

Oi 

CO 

tH 

eg 

CO 

4 

CO 

fH 

eg 

4 

4 

4 

4 

4 

CO 

eo 

QO 

flO 

rs. 

CD 

CD 

eg 

to 

fH 

40 

rg 

•rT 

*4* 

4* 

4 

CO 

CO 

CO 

CO 

eg 

eg 

eg 

fH 

eo 

eo 

(O 

eo 

eg 

eg 

eg 

eg 

fH 

rH 

rH 

O 

O 

O 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

Q 

O 

O 

O 

O 

< 

M 

M 

M 

w 

M 

w 

M 

M 

M 

w 

M 

M 

Vi 

u 

Ui 

Ui 

Ui 

Ui 

Ui 

Ui 

M 

IP 

IP 

IP 

IP 

M 

•4* 

O) 

eg 

fH 

O) 

40 

40 

CO 

00 

00 

CO 

Pg 

G) 

CD 

eg 

o 

Oi 

Oi 

4 

o 

o 

o 

o 

o 

X 

rx. 

O) 

CD 

C7> 

eg 

o 

rs. 

4 

CO 

CO 

o 

CD 

w 

rs. 

rs. 

eg 

eg 

4 

A 

4 

4 

o 

o 

o 

o 

u 

o 

O) 

o 

CO 

CO 

00 

CO 

G) 

a> 

40 

00 

Oi 

CO 

CO 

CO 

00 

eo 

G> 

o 

o 

o 

o 

CO 

O) 

40 

eg 

fH 

CO 

0) 

00 

fH 

40 

4 

CD 

CO 

cn 

o 

rs. 

CD 

eo 

cn 

4 

o 

o 

o 

o 

fx.. 

eg 

fH 

CD 

O 

00 

O) 

G> 

rg 

eg 

4 

rg 

4 

cn 

r^ 

a> 

fH 

eo 

cn 

CO 

r^ 

CO 

<4* 

CO 

CD 

fH 

eg 

4 

fH 

eg 

4 

fH 

eg 

eg 

eg 

eg 

eg 

rg 

fH 

fH 

iH 

rg 

eg 

eg 

eo 

eo 

CO 

CO 

4 

4 

in 

in 

CD 

CO 

o 

o 

o 

o 

o 

o 

O 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

O 

o 

o 

o 

o 

{.) 

M 

M 

M 

w 

M 

M 

u 

u 

u 

M 

M 

w 

2 

Ui 

M 

Ui 

Ui 

M 

M 

M 

Ui 

IP 

IP 

IP 

IP 

H-l 

O 

<• 

eg 

O 

o 

G) 

CO 

Oi 

eg 

40 

CD 

< 

rg 

eg 

QO 

4 

00 

4 

eg 

4 

eo 

eg 

CO 

eg 

2 

40 

(7> 

4- 

CD 

N 

CD 

O 

CO 

rs. 

r^ 

CD 

o 

> 

CD 

4 

rs. 

O 

OO 

rs. 

eo 

fs. 

Oi 

r< 

o 

4 

O 

•4* 

40 

CO 

CD 

O 

o 

o 

rs 

CO 

rs. 

40 

m 

O 

eo 

O 

eo 

4 

fs. 

OQ 

4 

Oi 

iH 

CO 

fs. 

rs. 

fH 

00 

CD 

CO 

eg 

CD 

ID 

tD 

eg 

Oi 

4 

eg 

4 

rs. 

rs. 

rs. 

a> 

O 

r^ 

fs. 

O) 

4 

H 

•4* 

o> 

es. 

O) 

O) 

40 

CO 

4 

O 

4 

eg 

u 

fH 

Oi 

to 

eo 

rg 

O) 

CO 

CD 

rH 

cn 

fH 

9 

CO 

CO 

4- 

40 

00 

eg 

eg 

eg 

CO 

CO 

4 

> 

o 

G 

4i 

o 

o 

o 

o 

N 

N 

eg 

eg 

• 

fH 

rg 

fH 

fH 

rg 

fH 

fH 

fH 

eg 

N 

eg 

0 

G 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

9 

a: 

O 

O 

o 

O 

O 

o 

O 

CO 

O 

o 

o 

o 

> 

•H 

9 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

eg 

< 

U 

1 

1 

9 

•0 

U 

> 

M 

M 

M 

w 

W 

w 

w 

w 

w 

M 

M 

Ui 

* 

> 

60 

> 

w 

w 

Ui 

M 

M 

Ui 

u 

IP 

IP 

IP 

IP 

IP 

a: 

•H 

O 

9 

< 

eg 

eg 

CO 

eg 

v\ 

CD 

4 

eg 

40 

CD 

CO 

oo 

* 

o 

CO 

< 

o 

o 

o 

O 

O 

o 

Oi 

r^ 

in 

Oi 

CO 

0 

•H 

9 

04 

cn 

eg 

fH 

rs. 

eg 

o 

rs. 

CO 

G> 

CO 

00 

4 

> 

B 

A 

B 

Q 

o 

o 

o 

O 

O 

o 

Oi 

CD 

Oi 

A 

eg 

eo 

«M 

o 

fH 

•4- 

O) 

G) 

O) 

eg 

rs. 

G) 

eg 

Oi 

CO 

9 

.H 

G 

X 

0 

o 

o 

o 

o 

o 

o 

Oi 

ID 

A 

fs. 

cn 

rH 

4^ 

o 

«M 

U 

H 

*H 

Oi 

sT 

O) 

4 

CO 

40 

G) 

G> 

tD 

CO 

<o 

A 

2: 

•H 

60 

•H 

o 

o 

o 

o 

o 

o 

Oi 

cn 

O 

o> 

cn 

O 

M 

CVJ 

c 

U 

O 

3 

CO 

O) 

•4- 

4 

eg 

eg 

CO 

CD 

oo 

CO 

CD 

CD 

G 

-0 

ig 

9 

fri 

CD 

> 

9 

0 

B 

> 

9 

Td 

9 

> 

9 

cn 

cn 

cn 

cn 

cn 

cn 

4 

eg 

rH 

4 

eg 

rH 

W 

9 

> 

R 

e 

<M 

U 

CO 

4’ 

40 

rs. 

eg 

4 

fs. 

fH 

eg 

4 

CD 

9 

p 

> 

A 

•0 

o 

G 

rH 

H 

ai: 

9 

<•> 

3 

O 

o 

s 

60 

O 

9 

G 

3 

rg 

9 

rg 

c/a 

a: 

9 

CD 

B 

•H 

rH 

rH 

eg 

eg 

eg 

eg 

eg 

CO 

4 

4 

CO 

40 

9 

a: 

9 

fH 

•H 

o 

fH 

eg 

eg 

eg 

eg 

eg 

eo 

4 

4 

cn 

cn 

2 

60 

9 

M 

O 

4) 

O 

O 

o 

o 

o 

o 

o 

o 

O 

O 

o 

O 

9 

60 

B 

O 

> 

G 

u 

O 

o 

o 

o 

o 

o 

o 

o 

O 

O 

o 

o 

9 

•H 

40 

1 

A 

3 

G 

9 

O 

1 

1 

4i 

C 

-G 

(M 

A 

eg 

4^ 

« 

01 

u 

u 

u 

u 

u 

u 

M 

w 

u 

Ui 

M 

Ui 

G 

M 

• 

G 

9 

•H 

> 

> 

■H 

9 

pg 

w 

Ui 

M 

Ui 

Ui 

M 

Ui 

Ui 

IP 

IP 

IP 

IP 

O 

O 

c 

9 

G 

0 

1 

9 

W 

40 

4' 

CD 

CD 

0) 

4 

rs. 

CO 

eg 

Oi 

CO 

CO 

G 

9 

o 

A 

9 

9 

9 

A 

A 

W 

cn 

CO 

CD 

4 

CD 

cn 

eg 

eg 

cn 

eo 

Oi 

< 

o 

o 

u 

9 

c 

3 

3 

o 

H 

rH 

40 

rs. 

40 

CO 

CD 

4 

eg 

eg 

4 

CO 

CD 

9 

R 

9 

G 

B 

B 

c a: 

9 

60 

JH 

CO 

eg 

cn 

Oi 

eg 

Oi 

eo 

Oi 

4 

fH 

cn 

fs. 

o 

u 

«2> 

u 

« 

4^ 

C 

B 

CO 

eo 

eg 

40 

O) 

eg 

o> 

CO 

G> 

4 

rs. 

CO 

rs. 

B 

60 

U 

9 

3 

G 

O 

A 

A 

B 

CO 

rs. 

fH 

rs. 

4 

o 

Oi 

fs. 

4 

O 

CD 

CO 

eg 

u-> 

«0 

9 

3 

u 

G 

n 

o 

B 

rv 

rs. 

4 

o 

o> 

rs. 

4 

O 

CD 

CD 

eg 

60 

9 

A 

R 

9 

•H 

•H 

•H 

•H 

eg 

iH 

fH 

CO 

o 

00 

o 

eg 

4 

OO 

m 

O) 

9 

n 

9 

3 

9 

4> 

0 

N 

OO 

o 

OO 

o 

eg 

4 

00 

rg 

CO 

9 

9 

60 

40 

A 

B 

Q 

CO 

» 

9 

9 

0 

60 

4^ 

B 

9 

<M 

9 

19 

O 

p 

9 

K 

CO 

Oi 

CD 

4 

fH 

lO 

CO 

r4 

CD 

o 

(M 

« 

a* 

« 

4d 

9 

»H 

fH 

O) 

to 

4 

fH 

fH 

CO 

to 

fH 

CO 

fH 

CM 

O 

10 

B 

IQ 

9 

9 

K 

9 

fl 

04 

M 

o 

<M 

•0 

M 

60 

M 

M 

0 

•ri 

0 

U 

A 

•H 

04 

o 

O 

« 

9 

cr 

>» 

O 

O 

o 

o 

O 

O 

O 

Oi 

Oi 

4 

rg 

O 

o 

60 

M 

U 

B 

40 

fH 

fH 

rg 

fH 

eg 

eg 

eg 

eg 

eo 

CO 

CO 

4 

m 

«o 

U 

Ed 

> n 

•o 

«0 

O 

4^ 

O 

o 

o 

o 

O 

O 

o 

o 

o 

40 

fs. 

rs. 

9 

60 

B 

Q 

9 

0 

• 

9 

S 

10 

M 

O 

O 

o 

O 

o 

o 

o 

o 

o 

O 

O 

O 

M 

c 

« 

9 

•H 

o 

o 

o 

o 

O 

o 

o 

40 

00 

0) 

fs. 

rs. 

U 

B 

•H 

9 

«a 

>% 

U 

A 

9 

M 

1 

•H 

c 

M 

M 

a> 

4^ 

> 

40 

40 

40 

40 

40 

40 

CO 

o 

00 

CD 

Oi 

G) 

•H 

B 

M 

0 

60 

or 

9 

B 

o 

U 

Ui 

M 

Ui 

w 

Ui 

Ui 

Ui 

M 

IP 

IP 

IP 

IP 

c 

•H 

4^ 

9 

1 • 

9 

•H 

< 

r*» 

rs. 

fs. 

rs. 

rs. 

rs. 

rs. 

o 

Oi 

G) 

4 

fH 

g 

B 

4^ 

fH 

U 

60  «Q 

s 

^4 

9 

O 

eg 

4 

rs. 

in 

rs. 

CD 

o 

CO 

eo 

CO 

eg 

fs. 

c 

0 

4J 

n 

9 

A 

H 

CM 

g 

fg 

9 

B 

3 

o 

A 

n 

m 

cn 

CD 

eo 

eg 

o 

CO 

cn 

eg 

CD 

fs. 

o 

O 

G*  <M 

U 

N O 

3 

3 

O) 

C7> 

O) 

O) 

O) 

a> 

G) 

40 

fH 

O 

O 

O 

o 

.H 

60  CM 

>s 

G 

<9 

40 

M 

A 

X 

9 

CD 

o 

CD 

eo 

to 

eo 

Oi 

oo 

eo 

rH 

CD 

eg 

«o  *o 

0 

9 

60 

60  « S 

4J 

•0 

« 

J3 

-4* 

O) 

4 

O) 

4 

60 

9 

O 

60 

9 

60 

9 

60 

i 

rs. 

eg 

cn 

00 

00 

fH 

o 

o 

eo 

iH 

o 

eo 

M 

9 

X 

9 

9 

« 9 

0 

B 

0 

eg 

u 

9 

A 

M 

o 

o 

04  -H 

M 

3 

9 

u 

*C 

o 

•o 

•o 

U 

tr 

£ 

0 

M 

9 

A 

A 

<9 

CM 

IQ 

1 

M 

9 

u 

eo 

eo 

eg 

A 

a> 

4 

eg 

fH 

4 

eg 

4 

«H 

N 

4^ 

a 

< a 

(0 

D 

B 

04 

A 

A 

B 

0 

•H 

fH 

60 

K 

G 

9 

0L4 

O) 

R 

60 

* 

M 04  « 

R 

A 

9 

60 

9 

O 

B 

9 

04 

W 

««• 

3 

« 

« 

c 

60 

m 

rH 

o 

O 

o 

o 

O 

o 

O 

eg 

fs. 

o 

fs. 

rs. 

3 

9 

B 

9 

A 

s..^ 

•H 

CM 

o 

o 

o 

O 

o 

O 

o 

eg 

rs. 

o 

fs. 

fs 

l-» 

c 

o 

G 

O 

« 

u o 

o 

O 

M 

« 

o 

O 

o 

o 

o 

o 

o 

4 

O 

4 

4 

rs. 

B 

eg 

9 

9 

O 

1 

CO 

fH 

0 

> 

o 

o 

o 

o 

o 

O 

o 

4 

o 

4 

4 

fs 

CO 

>% 

>v  ^ 

•H 

•G  .P 

M 

A 

CD 

9 

> 

o 

o 

o 

o 

o 

o 

o 

CO 

Oi 

4 

O 

eg 

>v  « 

rg 

60 

B 

M 

04 

60 

o 

o 

o 

o 

o 

o 

O 

o 

eo 

Oi 

4 

O 

eg 

2 

K 

00 

OC 

G4 

3 3 

U 'O 

3 

9 

•H 

l-t 

o 

o 

o 

o 

o 

o 

o 

CO 

CO 

Oi 

fH 

o 

K 

60  « 

9 

9 

u 

Ui 

60 

0 

fH 

1-4 

o 

o 

o 

o 

o 

O 

o 

eo 

CO 

Oi 

iH 

o 

9 

M 

M 

u* 

« 

« 

M -H 

<0 

> 

o 

o 

o 

o 

o 

o 

o 

fH 

o 

Oi 

CO 

eg 

9 

U 

9 

A 

M 

•H 

B 

* 

9 

•H 

•H 

•H 

o 

o 

o 

o 

o 

O 

o 

fH 

o 

Oi 

cn 

eg 

n 

9 

« 

4J 

0 

m ^ 

•H 

O 

M 

*0 

9 

A 

A 

CD 

u 

3 

CO 

U 

A 

K 

> 

Qli 

o 

c 

G 

C 

« 

4) 

u i2 

.c 

04 

3 

3 

o 

o 

o 

o 

o 

o 

o 

40 

eg 

o 

O 

o 

B 

G 

9 

9 

> 

9 

rg 

A 

9 

9 

CO 

o 

o 

o 

o 

o 

O 

o 

cn 

eg 

o 

o 

o 

U4 

m 

9 

• 

o< 

10 

u 

u 

u 

B 

10 

•o 

o 

40 

o 

40 

eg 

fH 

9 

A 

04 

C0 

> 

(Q 

o 

9 

U 

B 

> 

in 

o 

cn 

o 

cn 

eg 

A 

eg 

eg 

fH 

eg 

eg 

fH 

fH 

CO 

H 

o 

2 

H 

> 

(l4  > 

CJ 

< 

PQ 

W 

O 

> 

2 

H 

<7 

pg  > 

2 

CO 

eg 

M 

H 

< 

< 

w 

< 

»-l 

O 

<4 

w 

2 

2 

fH 

•H 

fH 

fH 

4 

eg 

CD 

Oi 

m 

Ui 

< 

< 

04 

H 

CO 

< 

pg 

2 

rg 

rg 

fH 

fH 

fH 

rg 

A 

4 

eg 

CD 

Oi 

fs 

H 

tH 

04 

m 

E° 

U 

< 

5 

r^ 

»H 

rs. 

G> 

O 

eg 

CO 

4 

CD 

H 

Q 

> 

Ui 

A 

X 

2 

5i 

iH 

fH 

fH 

fH 

r^ 

Oi 

O 

eg 

CO 

4 

CO 

CO 

o 

u 

0 

u 

< 

H 

CO 

eg 

CO 

4 

4 

4 

40 

CO 

40 

CO 

40 

CO 

Ui 

M 

rg 

N 

eo 

4 

4 

4 

in 

cn 

cn 

cn 

cn 

H 

CO 

C) 

H 

o 

ro 

CO 

CO 

< 

X 

28 


Table  3.  Sample  output  from  program  PTSUM 


Input  file:  PTRAN3D.TST 
Output  file:  PTSUM.TST 

PARTITION  OF  ENERGY  LOSS  IN  ENTIRE  PHANTOM: 

(a)  Initial  proton  energy:  160.000  MeV 

(b)  Energy  lost  In  Coulomb  Interactions:  lA 1.929  MeV 

Contribution  from  boundary  crossing  scores:!  1A1.813  MeV 
Contribution  to  from  track  ends ; 0 . 116  MeV 

(c)  Energy  lost  In  non-elastic  nuclear  Interactions;  18.06A  MeV 

(d)  Sum  of  (b)  and  (c):  159.993  MeV 

PARTITION  OF  ENERGY  LOSS  AT  VARIOUS  DEPTHS: 


z/r 


C( cross) 
C ( end ) 


C 

N 


depth  In  units  of  CSDA  range 
r - 1.765E+01  g/cin2 

Coulomb  Interactions,  crossing  of  scoring  boundaries 
Coulomb  Interactions,  track  ends 
C ( cross )+C( end) 

non-elastic  nuclear  Interactions 


Average  Energy  Loss  per  Unit  Depth,  MeV  cm2/g 


z/r 

CCcross ) 

C(end) 

C 

N 

C+N 

C(end)/C 

N/(C+N) 

0.000 

5.207 

0.000 

5.207 

1.592 

6.798 

0.00000 

0.23A13 

0.100 

5.333 

0.000 

5.333 

1.A83 

6.816 

0.00000 

0.21761 

0.200 

5.A91 

0.000 

5.  A91 

1.375 

6.866 

0.00000 

0.20026 

0.300 

5.692 

0.000 

5.692 

1.267 

6.959 

0.00000 

0.18205 

O.AOO 

5.952 

0.000 

5.952 

1.158 

7.111 

0.00000 

0.16291 

0.A50 

6.113 

0.000 

6.113 

l.lOA 

7.217 

0.00000 

0.15297 

0.500 

6.300 

0.000 

6.300 

1.0A9 

7.350 

0.00000 

0.1A276 

0.550 

6.522 

0.000 

6.522 

0.99A 

7.516 

0.00000 

0.13228 

0.600 

6.787 

0.000 

6.787 

0.939 

7.726 

0.00000 

0.12150 

0.650 

7.112 

0.000 

7.112 

0.882 

7.99A 

0.00000 

0.11036 

0.700 

7.520 

0.000 

7.520 

0.82A 

8.3AA 

0.00000 

0.09877 

0.720 

7.71A 

0.000 

7.71A 

0.800 

8.515 

0.00000 

0.09397 

0.7A0 

7.932 

0.000 

7.932 

0.776 

8.707 

0.00000 

0.08907 

0.760 

8.177 

0.000 

8.177 

0.750 

8.927 

0.00000 

0.08A06 

0.780 

8.A55 

0.000 

8.A55 

0.72A 

9.180 

0.00000 

0.07890 

0.800 

8.776 

0.000 

8.776 

0.697 

9.A73 

0.00000 

0.07360 

0.820 

9.150 

0.000 

9.150 

0.669 

9.819 

0.00000 

0.06811 

0.8A0 

9.593 

0.000 

9.593 

0.639 

10.232 

0.00000 

0.062A1 

0.860 

10.131 

0.000 

10.131 

0.606 

10.737 

0.00000 

0.056A5 

0.880 

10.801 

0.000 

10.801 

0.571 

11.372 

0.00000 

0.05020 

0.900 

11.668 

0.000 

11.668 

0.532 

12.201 

0.00000 

0.0A36A 

0.910 

12.210 

0.000 

12.210 

0.512 

12.722 

0.00000 

0.0A025 

0.920 

12.85A 

0.000 

12.85A 

0.A91 

13.3A5 

0.00000 

0.03679 

0.930 

13.638 

0.000 

13.638 

0.A68 

1A.107 

0.00000 

0.03320 

0.9A0 

1A.625 

0.000 

1A.625 

0.AA3 

15.068 

0.00000 

0.02939 

0.950 

15.930 

0.000 

15.930 

0.A12 

16.3A2 

0.00000 

0.02523 

0.955 

16.772 

0.000 

16.772 

0.39A 

17.166 

0.00000 

0.02295 

0.960 

17.795 

0.000 

17.795 

0.373 

18.168 

0.00002 

0.02051 

0 . 965 

19.117 

0.001 

19.119 

0.3A7 

19.A66 

0.00007 

0.0178A 

0.970 

20.875 

0.006 

20.881 

0.317 

21.198 

0.00027 

0.01A96 

0.975 

23.075 

0.018 

23.093 

0.281 

23.37A 

0.00078 

0.01202 

0.980 

25.802 

0.0A7 

25.8A9 

0.238 

26.086 

0.00181 

0.00912 

0.985 

28.353 

0.098 

28.A52 

0.189 

28.6A1 

0.003A5 

0.00660 

0.990 

29.5A9 

0.166 

29.71A 

0.138 

29.853 

0.00557 

0.00A63 

0.995 

27.866 

0.22A 

28.090 

0.091 

28.181 

0.00799 

0.00322 

1.000 

23.525 

0.2A5 

23.769 

0.052 

23.821 

0.01029 

0.00218 

1.005 

16.679 

0.21A 

16.893 

0.026 

16.919 

0.01268 

0.00151 

1.010 

9.919 

0.151 

10.070 

0.011 

10.081 

0.01A99 

O.OOlOA 

1.015 

A.87A 

0.085 

A.  960 

0.003 

A.  963 

0.01723 

0.00070 

1.020 

1.860 

0.039 

1.899 

0.001 

1.900 

0.020A9 

0.000A9 

1.025 

0.5A6 

O.OIA 

0.560 

0.000 

0.560 

0.025A3 

0.00036 

1.030 

0.130 

O.OOA 

0.13A 

0.000 

0.13A 

0.03123 

0.00022 

1.035 

0.016 

0.001 

0.017 

0.000 

0.017 

0.057A6 

0.00058 

l.OAO 

O.OOA 

0.000 

O.OOA 

0.000 

O.OOA 

0.0A319 

0.00000 
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Table  4. 


Fraction  of  energy  loss  (in  Coulomb  interactions)  due  to  track  ends,  as  fiinaion  of 
depth  in  units  of  CSDA  range,  zItq,  for  proton  beams  incident  with  energy  Tg. 


To,  MeV 
ro,  g/cm^ 

250.0 

37.94 

200.0 

25.96 

160.0 

17.65 

130.0 

12.28 

100.0 

7.718 

70.0 

4.080 

50.0 

2.227 

z/ro 

F R 

A C T I ( 

0 N 

0.950 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00001 

0.955 

0.00000 

0.00000 

0.00000 

0.00000 

0.00001 

0.00002 

0 . 00004 

0.960 

0.00001 

0.00001 

0.00002 

0.00003 

0.00004 

0.00009 

0.00018 

0.965 

0.00003 

0.00005 

0.00008 

0.00011 

0.00018 

0.00034 

0.00060 

0.970 

0.00013 

0.00019 

0.00028 

0.00039 

0.00058 

0.00100 

0.00166 

0.975 

0.00043 

0.00059 

0.00080 

0.00106 

0.00152 

0.00245 

0.00383 

0.980 

0.00108 

0.00142 

0.00186 

0.00238 

0.00325 

0.00499 

0.00745 

0.985 

0.00217 

0.00276 

0.00351 

0.00439 

0.00583 

0.00861 

0.01247 

0.990 

0.00362 

0.00452 

0.00563 

0.00693 

0.00903 

0.01301 

0.01854 

0.995 

0.00525 

0.00646 

0.00795 

0.00972 

0.01255 

0.01785 

0.02499 

1.000 

0.00690 

0.00844 

0.01034 

0.01255 

0.01607 

0.02263 

0.03147 

1.005 

0.00855 

0.01038 

0.01265 

0.01533 

0.01957 

0.02746 

0.03831 

1.010 

0.01018 

0.01228 

0.01502 

0.01817 

0.02307 

0.03232 

0 . 04446 

1.015 

0.01209 

0.01447 

0.01764 

0.02122 

0.02668 

0.03733 

0.05134 

1.020 

0.01430 

0.01694 

0.02043 

0.02462 

0.03133 

0.04300 

0.05947 

1.025 

0.01716 

0.02059 

0.02440 

0.02890 

0.03698 

0.05058 

0.06755 

1.030 

0.02197 

0.02632 

0.03148 

0.03551 

0.04471 

0.05934 

0.07954 

1.035 

0.03446 

0.03256 

0.04153 

0.04684 

0.05531 

0.07092 

0.09662 
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Table  5.  Energy  loss  from  nuclear  interactions,  as  fraction  of  the  energy  loss  from  Coulomb 
interactions,  as  function  of  depth  in  units  of  the  CSDA  range,  z/rQ,  for  proton  beams 
incident  with  energy  Tq. 


To.  MeV 
ro,  g/cm^ 

250.0 

37.94 

200.0 

25.96 

160.0 

17.65 

130.0 

12.28 

100.0 

7.718 

70.0 

4.080 

50.0 

2.227 

z/ro 

F R 

A C T I ( 

3 N 

0.000 

0.61257 

0.43368 

0.30570 

0.22122 

0.14880 

0.08914 

0.05524 

0.100 

0.55653 

0.39407 

0.27809 

0.20180 

0.13644 

0.08204 

0.05079 

0.200 

0.50000 

0.35411 

0.25037 

0.18233 

0.12398 

0.07477 

0.04623 

0.300 

0.44281 

0.31380 

0.22253 

0.16279 

0.11137 

0.06730 

0.04162 

0.400 

0.38489 

0.27320 

0.19458 

0.14316 

0.09852 

0.05957 

0.03695 

0.500 

0.32621 

0.23235 

0.16650 

0.12338 

0.08530 

0.05151 

0.03214 

0.600 

0.26689 

0.19125 

0.13827 

0.10319 

0.07154 

0.04313 

0.02696 

0.700 

0.20700 

0.14986 

0.10956 

0.08217 

0.05695 

0.03454 

0.02112 

0.800 

0.14653 

0.10773 

0.07941 

0.05959 

0.04123 

0.02504 

0.01433 

0.900 

0.08363 

0.06196 

0.04559 

0.03437 

0.02371 

0.01297 

0.00637 

0.950 

0.04739 

0.03514 

0.02583 

0.01881 

0.01182 

0.00546 

0.00207 

0.980 

0.01940 

0.01366 

0.00916 

0.00597 

0.00318 

0.00098 

0.00011 

0.990 

0.01068 

0.00725 

0.00463 

0.00286 

0.00137 

0.00031 

0.00002 

1.000 

0.00556 

0.00362 

0.00218 

0.00125 

0.00052 

0.00008 

0.00000 

1.010 

0.00299 

0.00184 

0.00104 

0.00055 

0.00018 

0.00002 

0.00000 

1.020 

0.00171 

0.00097 

0.00051 

0.00024 

0.00006 

0.00000 

0.00000 

1.030 

0.00097 

0.00059 

0.00023 

0.00010 

0.00003 

0.00000 

0.00000 

1.035 

0.00146 

0.00066 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 
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Table  6.  Comparison  of  the  average  energy  loss  (due  to  Coulomb  interactions)  per  unit  depth, 
(dE/^)c,  calculated  with  Grids  1,  2,  3,  4 and  5.  The  results  are  for  a 160-MeV 
beam  in  water,  and  for  each  case  are  based  on  a sample  of  1 million  Monte  Carlo 
histories. 


z/ro 

(dE/dz)c 
MeV  cm^/g 

Std.Dev. 

(Percent) 

Difference  w.r.t.  Grid 
(Percent) 

5 

(1) 

(2) 

(3) 

(4) 

Grid  1 

(5) 

Grid  2 

(6) 

Grid  3 

(7) 

Grid  4 

0.50 

6.300 

<0.01 

0.01 

0.01 

0.01 

<0.01 

0.90 

11.667 

<0.01 

0.06 

0.08 

0.05 

0.03 

0.95 

15.926 

0.01 

0.10 

0.14 

0.11 

0.07 

0.99 

29.459 

0.10 

0.47 

-0.18 

0.14 

-0.03 

1.00 

23.435 

0.16 

0.28 

-0.19 

-0.60 

-0.06 

1.01 

9.981 

0.30 

-0.97 

-1.29 

-2.02 

-0.71 

1.02 

1.891 

0.78 

-1.42 

-2.48 

-2.70 

-1.54 

Column  (1):  Depth  in  units  of  CSDA  range  rQ.  At  160  MeV,  rQ  = 17.65  g/cm^. 

Column  (2):  Average  energy  loss  per  unit  depth  due  to  Coulomb  interactions,  calculated  with 
Grid  5.  Track-end  contribution  is  omitted. 

Column  (3):  Percent  standard  deviation  of  results  in  column  (2). 

Columns  (4)  to  (7):  Percent  differences  of  results  calculated  with  Grids  1,  2,  3 and  4, 
respectively,  from  results  in  column  (2)  calculated  with  Grid  5. 
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Table  7.  List  of  files  stored  on  the  1.44-Mb  distribution  disk.  The  files  are  contained  in  four 
self-extracting  archive  files,  PARAM.EXE,  VPREP.EXE,  MPREP.EXE  and 
PTRAN.EXE. 


PARAM.EXE 

Self-extracting  archive  file;  expands  to  1005A08  bytes 

PARAMA.FOR 

Fortran  source  code 

COMPOS.  WAT 

Composition  data  for  water 

STOPRANG.WAT 

Stopping-power  and  remge  table  for  water 

OXFIT.DAT 

Total  nonelastic  nuclear  cross  section  for  oxygen 

PARAMA.TBl 

Output  table,  Grid  1 

PANAMA. TB2 

Output  table,  Grid  2 

PANAMA. TB3 

Output  table,  Grid  3 

PANAMA. TBA 

Output  table.  Grid  A 

PANAMA. TBS 

Output  table.  Grid  S 

PARAMA.AR2 

Output  array  for  VPREPA  and  MPREPA,  Grid  2 

PANAMA. AR3 

Output  array  for  VPREPA  and  MPREPA,  Grid  3 

PANAMA. ARA 

Output  array  for  VPREPA  and  MPREPA,  Grid  A 

PANAMA . PT2 

Output  array  for  PTRANID,  PTRANIDX  or  PTRAN3D,  Grid 

2 

PANAMA . PT3 

Output  array  for  PTRANID,  PTRANIDX  or  PTRAN3D,  Grid 

3 

PANAMA . PTA 

Output  array  for  PTRANID,  PTRANIDX  or  PTRAN3D,  Grid 

A 

VPREP.EXE 

Self-extracting  archive  file;  expands  to  3670A57  bytes 

VPREPA.FOR 

Fortran  source  code 

VCONA.FOR 

Fortran  source  code 

VSAMPA.FOR 

Fortran  source  code 

VSUMA.FOR 

Fortran  source  code 

VPREPA.AR2 

Output  file  for  VCONA,  Grid  2 

VPREPA.AR3 

Output  file  for  VCONA,  Grid  3 

VPREPA.AR2 

Output  file  for  VCONA,  Grid  A 

VCON4  produces  output  files  VPREPA.BR2,  VPREP4.BR3  or 
VPREP4.BR4  for  use  in  PTRAN3D  or  PTRANID. 


MPREP.EXE  Self-extracting  archive  file;  expands  to  29S3884  bytes. 


MPREP«.FOR 

MCONA.FOR 

MSAMPA.FOR 

THSET2 

MOLCl.COF 

M0LC2.C0F 

MFREFA.AR2 

MPREPA.AR3 

MPREP4.ARA 


Fortran  source  code 
Fortran  source  code 
Fortran  source  code 

Reduced  deflection  angles  for  Moliere  distribution 
Mollere  expansion  coefficients  f(l} 

Moliere  expansion  coefficients  f(2) 

Output  file  for  MCONA,  Grid  2 
Output  file  for  MCONA,  Grid  3 
Output  file  for  MCONA,  Grid  A 


MCONA  produces  output  files  MFREPA.BR2,  MPREPA.BR3  or 
MPREPA.BRA  for  use  in  PTRAN3D  or  PTRANID. 


PTRAN.EXE  Self-extracting  archive  file;  expands  to  310A28  bytes. 


PTRAN3D.F0R 

THSET2 

GAUSS.DAT 

PTRAN3D.TST 

XPTRAN3D . TST 

PTRANID. FOR 

PTRANIDX.FOR 

PTRANIDY.FOR 

PTSUM.FOR 

PTSUMTAB.TST 

PTSUMARR.TST 


Fortran  source  code 

Reduced  deflection  angles  for  Moliere  distribution 

Arrays  used  for  sampling  from  Gaussian  distribution 

Output  from  test  run 

Auxiliary  output  from  test  ran 

Fortran  source  code 

Fortran  source  code 

Fortran  source  code 

Fortran  source  code 

Output  table  from  test  run 

Output  arrays  from  test  run 
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Table  8.  Multiplicative  correction  factor  kjjp  used  in  Eq  (A1.8).  Values  in  parentheses  were 
obtained  by  interpolation.  This  t^le  corrects  misprints  in  Table  2.2  of  Berger  and 
Wang  (1988). 


Thomas -Fermi  Hartree-Fock  Potentials 

Potential 


Za/^ 

Z - 1 

2 

4 

6 

7 

8 

13 

29 

47 

0.00 

1.037 

0.863 

1.530 

1.129 

1.058 

1.081 

1.116 

1.186 

1.101 

1.183 

0.05 

(1.034) 

0.861 

1.527 

1.127 

(1.055) 

1.078 

1.114 

1.183 

1.098 

1.180 

0.1 

1.028 

0.855 

1.517 

1.119 

1.048 

1.071 

1.106 

1.175 

1.091 

1.172 

0.2 

1.004 

0.836 

1.482 

1.094 

1.024 

1.047 

1.081 

1.149 

1.066 

1.145 

0.4 

0.950 

0.790 

1.401 

1.034 

0.969 

0.990 

1.022 

1.086 

1.008 

1.083 

0.6 

0.918 

0.764 

1.355 

1.000 

0.937 

0.957 

0.989 

1.050 

0.975 

1.047 

0.8 

0.912 

0.759 

1.346 

0.994 

0.931 

0.951 

0.982 

1.043 

0.969 

1.040 

1.0 

0.918 

0.764 

1.356 

1.001 

0.937 

0.957 

0.989 

1.050 

0.975 

1.047 

1.2 

0.929 

0.773 

1.371 

1.012 

0.947 

0.968 

1.000 

1.062 

0.986 

1.059 

2.4 

0.968 

0.805 

1.428 

1.054 

0.986 

1.008 

1.041 

1.106 

1.027 

1.102 

4.8 

0.983 

0.818 

1.450 

1.070 

1.001 

1.023 

1.057 

1.122 

1.041 

1.117 

9.6 

0.987 

0.821 

1.454 

1.072 

1.003 

1.025 

1.058 

1.123 

1.040 

1.113 
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Figures 


Fig.  1.  Total  nonelastic  nuclear  cross  sections  for  protons  incident  on  oxygen- 16.  The  curve 
represents  the  cross  section  used  in  PTRAN.  The  experimental  points  are  from 
Renberg  et  al.  (1972)  and  Carlson  et  al.  (1975). 

Fig.  2.  Vavilov  energy-loss  distributions  for  protons  in  water.  Results  were  calculated  for 

initial  energies  and  path  lengths  from  Grid  4 (see  Section  2.2).  The  various  parameters 
have  the  following  values: 


Energy 

(MeV) 

Path  Length 
(g/cm^) 

(MeV) 

K 

25 

0.02282 

0.03780 

0.05123 

0.6878 

50 

0.04002 

0.03459 

0.09863 

0.3093 

100 

0.06850 

0.03184 

0.1834 

0.1388 

200 

0.1113 

0.02958 

0.3205 

0.06135 

Fig.  3.  Dependence  of  proton  energy  spectra  on  the  step-size  grid  used  in  PTRAN.  The 

bottom  panel  shows  the  track-length  distribution  differential  in  energy,  y(T,z),  from  a 
160-MeV  proton  beam  in  water,  at  a depth  z = 0.1  rQ,  where  rQ  is  the  CSDA  range 
(17.65  g/cm^).  This  histogram  was  calculated  with  step-size  Grid  5.  The  top  panel 
shows  the  percentage  amounts  by  which  values  of  y(T,z)  calculated  with  Grids  2,  3 
and  4 differ  from  the  results  obtained  with  Grid  5. 

Fig.  4.  Similar  to  Fig.  3,  for  a depth  z = 0.5  rQ. 

Fig.  5.  Similar  to  Fig.  3,  for  a depth  z = 0.9  rQ. 

Fig.  6.  Radial  distribution  of  energy  loss  due  to  Coulomb  interactions,  as  a function  of  the 

radial  distance  p from  a pencil  beam,  at  a depth  z = 0.1  rg,  for  a 160  MeV  beam  in 
water.  The  histograms  represent  Monte  Carlo  results  obtained  with  program 
PTRAN3D,  assuming  various  step-size  grids.  The  curve  is  calculated  using  the  theory 
of  Moii^re  (1955),  evaluated  in  the  continuous-slowing-down  approximation.  The 
Monte  Carlo  results  were  obtained  taking  into  account  multiple-scattering  angular 
deflections,  but  neglecting  lateral  multiple-scattering  displacements  within  each  step. 

Fig.  7.  Similar  to  Fig.  7,  but  calculated  with  inclusion  of  lateral  multiple-scattering 
displacement  in  each  step. 
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